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Abstract 

Current  methods  for  estimating  the  wavefront  slope  at  the  aperture  of  a  telescope  using  a 
Hartmann  wavefront  sensor  are  based  upon  a  centroid  shift  estimator.  The  centroid  shift  estimator 
determines  the  displacement,  or  shift,  of  the  centroid  off  the  optical  axis  using  a  moment  calculation 
of  the  intensity  distributions  recorded  in  each  subaperture.  This  centroid  shift  is  proportional  to  the 
average  slope  of  the  wavefront  in  each  subaperture.  A  maximum  a-posteriori  (MAP)  slope  estimator 
takes  advantage  of  a-priori  knowledge  of  the  wavefront  slope  statistics  and  total  irradiance  falling 
on  the  subaperture  detector  arrays  when  determining  the  shift  estimate.  In  order  to  derive  a  closed 
form  solution  for  the  MAP  estimator,  several  assumptions  were  made:  infinite  resolution  on  the 
detector  arrays,  no  read  noise  in  the  detection  process,  and  no  intensity  spillover  into  adjacent 
subapertures.  By  implementing  the  Hartmann  wavefront  sensor  and  MAP  estimator  in  simulation, 
the  performance  of  the  MAP  estimator  was  evaluated  using  realizable  wavefront  sensor  parameters. 
While  the  MAP  estimator  mean  square  error  (MSE)  performance  decreased  relative  to  the  centroid 
estimator  MSE  performance  as  a  result  of  spillover,  finite  detector  resolution,  and  read  noise,  the 
MAP  estimator  MSE  performance  was  found  to  be  upper  bounded  by  the  centroid  estimator  MSE 
in  all  cases. 
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EVALUATION  OF  A  MAXIMUM 
A-POSTERIORI  SLOPE  ESTIMATOR  FOR 
A  HARTMANN  WAVEFRONT  SENSOR 

1.  Background 

1.1  Introduction 

It  is  well  known  that  atmospheric  turbulence  reduces  the  resolution  of  imaging  systems.  In 
fact,  Isaac  Newton  recognized  that,  without  correction,  the  resolution  of  telescopes  was  limited 
not  by  aperture  size,  but  by  the  turbulence  in  the  atmosphere  [16].  A  variety  of  techniques  have 
been  used  to  improve  image  resolution.  One  solution  was  to  build  telescopes  at  higher  altitudes. 
While  the  higher  altitudes  significantly  reduced  the  effects  of  turbulence,  this  approach  had  its 
obvious  limitations.  Another  solution  was  to  compensate  for  the  turbulence  using  post-processing 
or  adaptive  optics  techniques.  While  it  has  been  shown  that  adaptive  optics  systems  can  improve 
resolution  by  compensating  for  the  atmospheric  effects,  cost  and  complexity  often  encourage  the 
use  of  speckle  imaging  [2, 29]  and  other  post-processing  techniques  such  as  inverse  filtering  or  blind 
deconvolution  [1,20].  Significant  research  has  been  done  by  the  Air  Force  Maui  Optical  Station 
(AMOS)  and  the  Starfire  Optics  Range  to  improve  upon  the  resolution  of  systems  imaging  through 
the  atmosphere  using  adaptive  optics  systems  [9,19].  One  of  the  key  components  of  the  adaptive 
optics  system  is  a  wavefront  sensor.  The  wavefront  sensor  gathers  information  used  to  estimate  the 
phase  of  the  incident  wavefront  so  the  phase  delays  that  distort  the  image  and  reduce  the  overall 
resolution  of  the  imaging  system  can  be  removed.  If  a-priori  information  about  the  atmospheric 
turbulence  is  used  along  with  the  information  gathered  by  the  wavefront  sensor,  better  estimates  of 
the  wavefront  phase  can  be  made.  This  thesis  examined  the  advantages  of  using  a-priori  knowledge 
about  the  atmospheric  turbulence  when  gathering  data  to  reconstruct  and  remove  the  incident 
wavefront  phase  distortions.  The  following  sections  present  the  information  required  to  understand 
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the  operation  of  a  Hartmann  wavefront  sensor  and  its  role  in  the  adaptive  optics  system,  the  purpose 
of  the  research,  summary  of  results,  and  a  brief  thesis  overview. 

1.2  Atmospheric  turbulence 

Atmospheric  turbulence  is  caused  by  the  heating  and  cooling  of  the  Earth’s  surface  by  the 
sun.  Uneven  temperature  distributions  in  the  atmosphere  result  in  the  formation  of  eddies  with 
varying  indexes  of  refraction.  These  random  refractive-index  inhomogeneities  distort  the  light  that 
propagates  any  significant  distance  through  the  atmosphere  [8].  Consider  light  propagating  from  a 
distant  star.  Initially,  the  light  propagates  outward  with  a  spherical  wavefront.  After  a  significant 
distance,  the  wavefront  can  be  modeled  as  planar.  As  the  plane  wave  propagates  through  the 
atmosphere,  different  parts  of  the  wavefront  pass  through  eddies  with  differing  refractive  indexes. 
As  a  result,  the  surface  of  the  wavefront  experiences  various  phase  delays  and  appears  to  be  dimpled. 
This  dimpling  effect  is  a  direct  result  of  the  light  passing  through  pockets  of  random  index  of 
refraction.  The  most  significant  effect  of  atmospheric  turbulence  is  that  it  imparts  a  random  tilt  on 
the  wavefront  [8].  For  long  exposure  images,  the  random  tilt  will  cause  the  image  to  roam  around 
the  image  plane  of  the  detector,  resulting  in  a  blurred  image. 

Atmospheric  turbulence  is  characterized  by  a  set  of  parameters  which  include  the  outer  scale, 
Lo,  and  the  Pried  parameter.  To-  The  outer  scale  is  a  measure  of  the  largest  turbulent  eddies  that 
are  of  a  single  index  of  refraction  [22].  The  Pried  parameter  is  the  atmospheric  coherence  diameter 
which  relates  to  the  overall  strength  of  the  turbulence  induced  perturbations  [11].  The  turbulent 
nature  of  the  atmosphere  is  modeled  using  wave  structure  functions.  A  wave  structure  function 
describes  the  phase  and  amplitude  fluctuations  of  a  wave  due  to  turbulence  [5].  For  the  purpose 
of  the  simulations  in  this  thesis,  phase  structure  function  statistics  will  only  be  considered  since 
the  effects  of  amplitude  perturbations  are  small  relative  to  the  effects  of  phase  perturbations  [3]. 
Two  common  phase  structure  functions  that  are  used  to  model  the  effects  of  turbulence  are  the 
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Figure  1.1  Wave  front  distortions  caused  by  atmospheric  turbulence 


Kolmogorov  structure  function  and  the  von  Karman  structure  function.  The  Kolmogorov  structure 
function  assumes  infinite  outer  scale  and  is  the  most  commonly  used  model.  The  von  Karman 
structure  function  uses  a  finite  outer  scale.  The  outer  scale,  Lo,  is  the  size  of  the  largest  turbulent 
eddie  that  is  of  a  single  index  of  refraction.  Turbulent  eddies  larger  than  Lg  are  not  assumed  to  be 
of  a  single  index  of  refraction,  but  composed  of  smaller  eddies  with  various  indexes  of  refraction. 

The  affect  atmospheric  turbulence  has  on  the  imaging  process  is  that  it  can  only  degrade  the 
resolution  of  the  imaging  system.  As  the  Fried  parameter,  Vg,  becomes  small,  the  turbulent  nature 
of  the  atmosphere  becomes  the  limiting  factor  in  the  angular  resolution  of  imaging  systems,  not 
the  aperture  size.  Fried  showed  that  when  the  telescope  diameter,  D,  is  much  less  than  Vg,  the 
diameter  is  the  limiting  factor.  But,  when  D  becomes  larger  than  rg,  the  resolution  is  consistent 
with  a  telescope  with  diameter  r<,  [7]. 
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Figure  1.2  A  typical  wavefront  compensation  adaptive  optical  system 
1.3  Adaptive  optics  systems 

The  fundamental  goal  of  adaptive  optics  systems  is  to  remove  the  perturbations  caused  by 
atmospheric  turbulence  in  real  time,  thus  improve  the  overall  resolution  of  the  imaging  system.  The 
removal  of  perturbations  is  accomplished  by  performing  two  basic  functions:  wavefront  sensing  and 
wavefront  correction.  A  wavefront  sensor  measures  the  turbulence  induced  phase  distortions  across 
the  telescope  aperture.  This  data  is  then  used  to  control  a  wavefront  modifying  device.  The  adaptive 
optics  system  shown  in  Fig.  1.2  uses  a  technique  known  as  wavefront  compensation  and  is  used  when 
the  objective  is  to  obtain  the  best  possible  image  of  a  distant  object  viewed  through  a  turbulent 
path  [13].  The  wavefront  sensor  provides  the  information  required  to  drive  a  deformable  mirror. 
The  purpose  of  the  deformable  mirror  is  to  remove  the  phase  perturbations  in  the  wavefront.  A 
tip-tilt  mirror  can  also  be  used  to  remove  the  overall  tilt  of  the  incident  wavefront  [6].  This  process 
is  known  as  global  tilt  removal.  Ideally,  the  wavefront  correction  would  completely  negate  the 
phase  perturbations  that  are  present  in  the  telescope  aperture.  However,  there  are  many  factors 
that  contribute  to  error  in  the  adaptive  optics  process  which  result  from  the  inability  to  produce 
perfect  wavefront  sensors  and  deformable  mirrors  [10]. 
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1.4  Hartmann  wave-front  sensor 


While  numerous  wavefront  sensing  techniques  exist  [13],  this  research  studies  the  performance 
of  a  Hartmann  wavefront  sensor.  The  Hartmann  wavefront  sensor  is  a  tilt  or  slope  sensor.  It  consists 
of  an  array  of  lenslets,  or  subapertures,  that  segments  the  incident  wavefront  in  the  pupil  of  the 
imaging  system.  An  intensity  distribution,  or  centroid,  is  formed  behind  each  subaperture  on  a 
charged  coupled  device  (CCD)  detector  array.  The  distance  the  centroid  is  displaced  from  the 
optical  axis  is  called  a  shift.  The  shift  is  directly  related  to  the  average  slope  of  the  wavefront  in 
the  pupil  of  the  subaperture  [3].  Using  the  centroid  shift  data  from  all  of  the  subapertures  in  the 
wavefront  sensing  array,  the  phase  of  the  incident  wavefront  can  be  reconstructed  [6,25,28].  The 
accuracy  of  the  shift  data  will  directly  affect  the  accuracy  of  the  wavefront  reconstruction  as  well 
as  the  ability  of  the  adaptive  optics  system  to  remove  the  phase  perturbations. 

1.5  Maximum  a-posteriori  estimators 

The  accuracy  of  the  data  collected  by  the  Hartmann  wavefront  sensor  is  limited  by  the  random 
nature  of  the  atmospheric  turbulence  and  the  photo  detection  process.  If  the  random  nature  of  the 
atmospheric  turbulence  and  the  detection  process  can  be  characterized  using  probability  density 
functions  (PDFs),  shift  estimators  that  incorporate  the  PDF  information  can  be  developed.  If 
the  wavefront  turbulence  statistics  are  known  a-priori,  the  atmospheric  turbulence  statistics  can 
be  characterized  using  a  posterior  PDF  and  a  maximum  a-posteriori  (MAP)  estimator  can  be 
developed.  The  approach  taken  by  Sallberg  [23]  was  to  model  the  image  detection  process  as 
a  product  of  joint  PDFs.  By  maximizing  the  joint  PDF  with  respect  to  the  shift  parameter,  a 
MAP  shift  estimator  was  derived.  Several  assumptions  were  made  in  order  to  obtain  a  closed  form 
solution  for  the  MAP  estimator: 

•  The  intensity  distribution  on  the  image  plane  was  modeled  as  a  Gaussian  distribution. 
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•  The  intensity  distribution  was  sufficiently  compact  and  the  local  wavefront  tilts  were  suffi¬ 
ciently  small  so  that  the  intensity  distribution  formed  by  one  subaperture  does  not  bleed  over 
onto  the  detector  pixels  from  an  adjacent  subaperture. 

•  The  CCD  detector  array  was  assumed  to  have  infinite  spatial  resolution. 

•  The  random  nature  of  optical  detection  was  modeled  with  Poisson  statistics. 

•  The  detection  process  had  no  read  noise. 

The  resulting  MAP  estimator  incorporated  statistical  knowledge  of  wavefront  slopes  and  the  de¬ 
tected  light  levels  when  estimating  the  shift  data.  If  all  the  parameters  in  the  joint  PDF  were 
jointly  Gaussian,  the  resulting  estimator  would  be  the  Minimum  Mean  Square  Error  (MMSE)  es¬ 
timator.  Since  the  arrival  of  photons  on  the  CCD  detector  was  modeled  as  a  Poisson  process  and 
the  intensity  distribution  was  only  approximated  by  a  Gaussian  distribution,  the  MAP  estimator 
will  not  be  an  MMSE  estimator.  Sallberg  found  that  the  MAP  estimator  mean  square  error  (MSE) 
had  the  centroid  estimator  MSE  as  its  upper  bound. 

1.6  Goal  of  research 

The  purpose  of  this  research  was  to  evaluate  the  performance  of  Sallberg’s  MAP  estimator  [23] 
under  actual  operating  conditions.  Many  of  the  assumptions  made  by  Sallberg  no  longer  hold 
when  considering  a  realizable  wavefront  sensor  design.  The  MAP  estimator  must  be  evaluated 
using  simulation  since  a  closed  form  solution  can  not  be  found  when  realizable  wavefront  sensor 
parameters  are  considered.  The  following  MATLAB  5.0©  simulations  were  used  to  complete  the 
research: 

•  A  Fourier  series  based  phase  screen  generator  modified  from  code  developed  by  Welsh  [30]. 

•  A  Hartmann  wavefront  sensor  employing  an  AT  x  iV  array  of  square  subapertures. 
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•  A  slope  correlation  matrix  generator  for  von  Karman  and  Kolmogorov  atmospheric  statistics 
for  use  in  the  MAP  estimator. 

1. 7  Overview 

Chapter  II  covers  the  principles  behind  Hartmann  wavefront  sensors  and  the  simulation  de¬ 
velopment.  Chapter  III  discusses  the  implementation  of  Sallberg’s  MAP  estimator  and  the  required 
equations  for  the  slope  correlation  matrix.  Chapter  IV  verifies  the  simulation  accuracy,  discusses 
the  various  parameters  used  in  the  test  simulations,  and  presents  a  test  matrix.  Chapter  V  exam¬ 
ines  the  performance  of  the  MAP  estimator  in  simulation  and  finally.  Chapter  VI  summarizes  the 
thesis. 
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II.  Imaging  theory  and  wavefront  sensor  simulation  development 

2.1  Introduction 

The  purpose  of  this  research  was  to  develop  a  computer  simulation  for  a  Hartmann  wavefront 
sensor  and  use  it  to  evaluate  the  performance  of  a  maximum  a-posteriori  (MAP)  slope  estimator. 
While  the  previous  chapter  presented  a  basic  overview  of  the  area  of  interest,  the  material  in  this 
chapter  is  focused  on  the  wavefront  sensor.  A  theoretical  foundation  for  using  a  linear  systems 
approach  to  model  the  imaging  process  of  a  single  lens  system  and  the  operation  of  the  Hartmann 
wavefront  sensor  are  explained.  Then,  the  computer  simulation  of  the  Hartmann  wavefront  sensor 
will  be  discussed  in  detail. 

2.2  Image  formation 

Light  propagating  in  a  source  free  media  system  can  be  analyzed  as  a  linear,  space-invariant 
(LSI)  system  using  Fourier  optics  theory  [12].  Using  the  LSI  framework,  the  single  lens  imaging 
system  with  an  image  plane  located  one  focal  length  away  pictured  in  Fig.  2.1  will  produce  an 
intensity  distribution  that  is  the  magnitude  squared  of  the  Fourier  transform  of  the  wavefront  in 
the  aperture. 


Wavefront 


Figure  2.1  Fourier  optics  approach  to  imaging  with  a  single  lens. 

In  an  actual  imaging  system,  both  the  wavefront  in  the  pupil  and  the  intensity  distribution 
on  the  detector  are  continuous.  In  order  to  implement  the  imaging  process  in  simulation,  the 


2-1 


wavefront  must  be  represented  using  a  finite  number  of  samples.  Representing  a  continuous  function 
with  a  finite  number  of  samples  makes  the  imaging  simulation  a  discrete  process,  and  a  Fast 
Fourier  Transform  (FFT)  of  the  wavefront  in  the  pupil  must  be  used  in  order  to  determine  the 
intensity  distribution  on  the  image  plane.  Two  problems  must  be  considered  when  working  in  the 
discrete  domain:  wrap  around  error  and  aliasing  due  to  under-sampling  [21].  In  order  to  avoid 
problems  associated  with  wrap  around  error  when  taking  FFTs,  the  wavefronts  are  placed  in  a 
zero-padded  array  before  transforming.  The  problem  of  under-sampling  is  handled  by  ensuring 
the  wavefront  sampling  is  larger-than  or  equal-to  the  sampling  that  is  used  by  the  phase  screen 
generation  program  [30]. 

2.3  Image  detection 

Each  subaperture  in  the  wavefront  sensor  uses  an  array  of  charged  coupled  device  (CCD) 
pixels  to  detect  the  intensity  distribution  that  forms  on  the  image  plane.  The  CCD  array  detects 
the  intensity  distribution  by  counting  the  arrival  of  photons  at  each  pixel  location.  Since  photons 
arrive  at  random  intervals,  the  detection  process  can  not  be  treated  as  deterministic.  In  addition  to 
a  random  arrive  rate,  the  detection  process  is  subject  to  several  sources  of  random  noise  including 
background  noise  and  readout  noise  [26]. 

2.3.1  Shot  noise.  Shot  noise  is  another  term  for  photo-conversion  noise.  It  accounts 
for  the  uncertainty  in  the  detection  of  photons.  Photons  arrive  at  a  random  rate  and  at  random 
locations  on  the  detector  array  which  was  modeled  by  a  Poisson  process.  A  Poisson  process  has 
a  variance  equal  to  its  mean  making  it  a  signal  dependent  process.  Since  shot  noise  is  signal 
dependent,  it  is  the  limiting  factor  in  the  detection  process  [28].  The  best  centroid  estimation  in 
the  presence  of  shot  noise  is  referred  to  as  shot  limited  performance  and  is  the  Cramer-Rao  lower 
bound  for  any  unbiased  shift  estimator  [32], 
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Figure  2.2  Geometry  for  a  single  subaperture  of  a  Hartmann  wavefront  sensor. 

2.3.2  Read  noise.  While  the  detection  of  the  photons  by  the  CCD  pixels  is  modeled  with 
a  signal  dependent  random  process,  the  act  of  reading  this  data  contributes  to  the  uncertainty 
by  adding  more  random  noise.  A  charge  coupled  device  array  operates  by  counting  the  arrival  of 
photons  in  each  pixel  over  a  given  period  of  time.  At  the  end  of  each  time  period,  the  photon  count 
is  read  out  of  the  detector  array.  The  random  noise  resulting  from  the  read  out  process  is  referred 
to  as  readout  noise  and  is  modeled  with  an  additive  zero  mean  Gaussian  random  process  [26]. 

2.4  Hartmann  wavefront  sensor 

The  Hartmann  wavefront  sensor  consists  of  an  array  of  lenslets  or  subapertures.  The  incident 
wavefront  in  the  aperture  of  the  imaging  system  is  segmented  over  the  array  of  subapertures.  The 
lenslet  in  each  subaperture  forms  an  intensity  distribution  onto  an  array  of  CCD  pixels  located  at 
the  lenslet  focal  length.  The  intensity  distribution  formed  by  each  lenslet  is  detected  by  an  array  of 
CCD  pixels,  and  the  centroid  shift  can  be  estimated  for  each  subaperture.  The  estimated  centroid 
shift  data  can  be  used  to  calculate  the  slope  of  the  wavefront  in  the  aperture  of  the  imaging  system. 
Referring  to  the  diagram  in  Fig.  2.2,  the  slope  of  the  wavefront  in  the  ith  subaperture,  can  be 
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Small  To 


Figure  2.3  The  wavefront  in  the  aperture  (a)  and  its  piecewise  linear  approximation  (b)  across 
an  array  of  subapertures  (c)  for  different  values  of  the  Pried  parameter,  To- 


determined  using  two  relationships.  First,  the  slope  is  proportional  to  the  centroid  shift,  Xi  [3]: 


radians 

meter 


\ 


(2.1) 


where  k  =  2ir/X,  X  is  the  average  wavelength,  and  /j  is  the  subaperture  focal  length.  Secondly, 
the  slope  of  the  wavefront  in  the  subaperture  can  be  found  using  a  relationship  developed  by 
Wallner  [28], 

=  (^^).  (2.2) 


where  VWi(f)  is  the  gradient  of  the  subaperture  weighting  function  of  the  ith  subaperture,  and 
<j){x)  is  the  wavefront  phase  in  the  ith  subaperture.  It  should  be  noted  that  the  expressions  defined 
in  Eqn.  (2.1)  approximates  the  wavefront  slope  as  linear  in  the  subaperture.  Since  the  Hartmann 
wavefront  sensor  is  an  array  of  subapertures,  the  wavefront  in  the  aperture  of  the  imaging  system 
is  averaged  in  a  piecewise  linear  fashion  [14].  When  the  Pried  parameter,  Tq,  is  large,  the  linear 
approximation  holds.  However,  as  Tq  decreases,  the  linear  approximation  will  not  be  valid  as  the 
spatially  averaged  wavefront  tilt  can  differ  greatly  from  the  true  wavefront  tilt  as  shown  in  Fig.  2.3. 
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2.5  Hartmann  wavefront  sensor  simulation 


In  order  to  develop  a  simulation  for  the  Hartmann  wavefront  sensor,  the  operation  must 
be  broken  up  into  several  steps.  The  approach  taken  to  model  the  operation  of  the  Hartmann 
wavefront  sensor  was  similar  to  one  described  by  Cannon  [3]  in  which  the  incident  wavefront  was 
segmented  over  an  array  of  microlenses  and  FFTs  were  used  to  determine  the  corresponding  inten¬ 
sity  distributions.  A  flow  chart  for  the  Hartmann  wavefront  sensor  simulation  is  shown  in  Fig.  2.4. 
The  simulation  input  is  a  phase  screen,  and  the  outputs  are  vectors  containing  the  actual  centroid 
shifts  (Xact)  and  the  wavefront  sensor  shift  estimates  The  actual  centroid  shift  can  be  calcu¬ 

lated  in  the  simulation  only  because  the  actual  phase  data  is  available.  An  actual  wavefront  sensor 
can  only  estimate  the  centroid  shift  based  upon  the  detected  signal  values.  The  difference  between 
the  actual  centroid  shift  and  the  estimated  centroid  shift  is  the  estimation  error.  The  simulation 
segments  the  incident  wavefront  over  an  array  of  subapertures,  calculates  the  actual  centroid  shifts, 
calculates  the  intensity  distribution  on  the  detector  array,  accounts  for  shot  noise,  read  noise,  and 
overlap  onto  adjacent  detectors,  then  estimates  the  centroid  shifts. 

The  major  steps  of  the  wavefront  sensor  simulation  are  described  in  detail  below: 

•  The  incident  wavefront  is  of  uniform  amplitude  and  varying  phase.  The  phase  of  the  wavefront  is 
the  input  to  the  simulation.  The  phase  data,  or  phase  screen,  of  the  incident  wavefront  is  segmented 
over  the  wavefront  sensor  subaperture  array.  The  actual  centroid  shift  is  then  calculated  for  each 
subaperture.  Using  Eqns.  (2.1)  and  (2.2),  the  actual  centroid  shift,  Sacu  is: 

Sact  =  -|  /  d^Wiix)cl>ix),  (2.3) 

where  /j  is  the  subaperture  focal  length,  Wi{x)  is  the  subaperture  weighting  function  of  the  ith 
subaperture,  k  =  27r/A,  and  <j>(x)  is  the  wavefront  phase  in  the  ith  subaperture.  Since  the  wavefront 
sensor  subapertures  are  assumed  square  in  the  simulation,  the  subaperture  weighting  function. 
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Figure  2.4  Flow  chart  for  the  Hartmann  wavefront  sensor  simulation 


Wi{x,y),  for  all  subapertures  can  be  modeled  as: 


Wiix,y) 


^rect{i)rect(i). 


(2.4) 


where  dw  is  the  width  of  the  wavefront  sensor  subaperture,  and 


rect 


1  <x<^ 

0  otherwise. 


(2.5) 
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Using  the  definition  of  the  gradient  operator,  V  =  ■§^x  +  the  gradient  of  the  subaperture 
weighting  function  is 


VWi{x,y)  = 


^Wi{x,y)x  +  ^Wiix,y)y 


(*(*'  t)  - 


(2.6) 


Now,  Eqn.  (2.3)  can  be  written  in  an  expanded  form  by  substituting  Eqn.  (2.6)  for  the  gradient  of 
the  subaperture  weighting  function: 


Xact  — 


+ 


;  f  y  d*d!,[(i(*  +  ^)  -  «(«  -  ^))  recl(A)i 

(«(„  +  -  ((y  -  rert(|.)s] . 


(2.7) 


An  expression  for  the  actual  centroid  shift  can  be  found  by  using  the  sifting  property  of  Dirac  delta 
functions  [12]: 


^aci  ~ 


^[/d!,(#(^.») *+/'i*('‘(“^'f)  -«'(». -f))  4  (2*) 


Consider  only  one  of  the  two  orthogonal  directions.  The  Xact  value  in  the  x  direction  can  be  written 
as 


^act 


A. 


X. 


(2.9) 
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Integrating  the  phase  with  respect  to  y  and  dividing  by  the  width  of  the  subaperture  results  in  the 
average  phase  value  at  the  points  x  ±  The  difference  between  the  average  phase  at  these  two 
points  scaled  by  fifkdw  is  the  actual  shift  in  the  x  direction.  In  the  discrete  domain,  the  actual 
shift  can  be  found  in  a  similar  manner.  In  order  to  determine  the  average  phase  value  at  the  points 
x±^,  the  phase  samples  are  summed  in  the  y  direction  at  the  points  x±^  and  divided  by  the 
total  number  of  samples.  The  difference  between  the  average  phase  at  these  two  points  scaled  by 
/( /kdw  is  the  actual  shift  in  the  x  direction: 


^act  — 


_j,  T.UK-^-y) 


kd^ 


N 


N 


X. 


(2.10) 


The  discrete  form  of  Eqn.  (2.8)  can  be  written  as 


^act  “ 


kd^[  N  N  ^ 


^  iV  N  ^ 


(2.11) 


where  is  the  total  number  of  samples  in  the  phase  screen,  //  is  the  subaperture  focal  length, 
k  =  2-k/X,  and  d^  is  the  width  of  the  wavefront  sensor  subaperture. 

•  After  the  actual  centroid  shifts  have  been  calculated  for  each  subaperture,  the  estimated 
centroid  shifts  are  determined.  The  first  step  is  to  zero-pad  each  segmented  piece  of  the  phase  screen 
in  order  to  minimize  the  effects  of  wrap  around  when  the  FFT  is  taken  [21].  The  zero-padding  is 
accomplished  by  resizing  the  segmented  piece  of  the  phase  screen  to  15  x  15  samples.  Based  upon 
the  selected  atmospheric  parameters,  the  sampling  of  the  phase  screens  generated  by  the  phase 
screen  generator  [30]  is  less  than  15  x  15  samples.  Resizing  to  a  larger  number  of  samples  ensures 
no  loss  of  data  or  aliasing  due  to  under-sampling.  The  phase  screen  data  is  placed  into  a  15  x  15 
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array  centered  on  a  larger  64  x  64  array  containing  0  values.  A  64  x  64  array  is  used  to  take 
advantage  of  FFT  algorithms. 

•  At  this  point  in  the  simulation,  the  phase  screen  in  each  subaperture  has  been  placed  into 
a  larger  zero-padded  array.  The  next  step  is  to  take  the  FFT  of  each  zero-padded  array.  The 
resulting  arrays  contain  complex  amplitude  data  for  each  subaperture.  As  shown  in  Fig.  2.5,  all  of 
the  complex  amplitude  data  is  added  to  a  single  detection  array.  Since  the  simulation  was  designed 
such  that  the  physical  size  of  the  arrays  containing  the  complex  amplitude  data  was  larger  than  the 
physical  size  of  the  subaperture,  the  complex  amplitude  data  for  adjacent  subapertures  will  overlap 
when  added  to  the  detection  array.  However,  if  overlap  between  subapertures  is  not  allowed,  an 
overlap  mask  is  placed  over  the  array  containing  the  complex  amplitude  data,  nulling  all  amplitude 
data  in  the  overlap  region.  Assigning  all  of  the  amplitude  values  in  the  overlap  region  to  zero  will 
negate  the  effect  of  overlap  on  the  detection  array. 

•  After  all  of  the  complex  amplitudes  are  added  to  the  detector  array,  the  detector  array  is 
multiplied  by  its  complex  conjugate  to  determine  the  intensity  distribution  on  the  detector  plane. 

•  Figure  2.6  shows  that  the  detection  array  contains  the  intensity  distributions  formed  by  all 
of  the  subapertures.  At  this  point  in  the  simulation,  the  detection  of  the  intensity  distribution  by 
the  CCD  array  is  modeled.  In  an  actual  wavefront  sensor,  each  subaperture  contains  an  N  x  N 
array  of  CCD  detector  pixels  located  in  the  image  plane.  In  the  simulation,  N  is  a  user  defined 
parameter.  The  CCD  array  is  modeled  as  a  square  array  of  pixels  with  no  spacing  between  them. 
To  account  for  the  effects  dead  space  between  the  detector  pixels  in  an  actual  CCD  array,  the 
photon  count  can  be  scaled  by  a  user  defined  duty  cycle  parameter.  The  intensity  distribution 
corresponding  to  each  subaperture  is  sampled  into  an  N  x  N  CCD  array.  The  sum  of  the  values  in 
each  subaperture  CCD  array  is  normalized  to  equal  the  average  photon  count  per  subaperture,  K. 
The  CCD  array  for  each  subaperture  now  contains  the  normalized  detected  signal  values. 
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Array  of  subapertures 


Figure  2.5  The  wavefront  is  segmented  over  an  array  of  subapertures.  The  individual  segments 
are  placed  in  zero-padded  arrays.  The  FFT  is  taken  to  get  the  complex  amplitude 
data,  an  overlap  mask  is  placed  on  the  data,  then  the  data  is  added  to  a  detection 
array. 


Figure  2.6  The  intensity  distribution  for  each  subaperture  is  sampled  into  a  smaller  CCD  array. 
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•  The  detection  process  by  a  CCD  is  subject  to  shot  and  read  noise.  To  account  for  the 
uncertainty  of  the  detection  process,  each  value  in  the  CCD  arrays  is  subjected  to  a  Poisson  process. 
The  CCD  arrays  now  contain  shot  limited  detected  values.  The  shot  limited  detected  values  in 
the  CCD  arrays  are  then  subjected  to  an  additive  Gaussian  zero  mean  process  to  account  for  the 
effects  of  read  noise.  The  variance  of  the  read  noise  is  another  user  defined  parameter. 

•  Now  that  the  CCD  arrays  contain  the  detected  signal  values,  the  estimated  centroid  shift 
can  be  calculated.  For  a  detector  array  with  a  finite  number  of  pixels,  the  centroid  shift  estimates 
for  each  subaperture,  jfi,  are  calculated  using  a  first  moment  calculation  [4], 


EiEjPij 


(2.12) 


where  Pij  is  the  photon  count  in  the  {i,j)  pixel  of  the  subaperture  detection  array,  and  Xi  and  yj 
are  the  offsets  from  the  center  of  the  subaperture  to  the  center  of  the  ith  or  jth.  pixel. 


2.6  Using  the  wavefront  sensor  simulation 

The  wavefront  sensor  simulation  is  set  up  as  a  MATLAB  5.0©  function  call.  The  code  listing 
in  Appendix  D  is  called  using  the  following  command: 


function  [xest,xact]  =  wfsl(PhaseScreen,  SubapertureArraySize,  SubapertureFocalLength, 
SubapertureDiameter,  AverageWavelength,  DetectorArraySize,  ReadNoise Variance,  ShotNoise, 
AveragePhotonCount,  WavefrontTilt,  DetectorDutyCycle,  SubapertureOverlap) 

The  wavefront  sensor  simulation  was  designed  to  allow  the  following  parameters  to  be  varied: 
•  PhaseScreen  -  The  phase  data  of  the  incident  wavefront  in  the  aperture  of  the  imaging  system. 
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•  SubapertureArraySize  -  The  squaxe  root  of  the  number  of  subapertures  in  the  wave  front 
sensor  array.  The  simulation  assumes  a  square  array  of  subapertures. 

•  SubapertureFocalLength  -  The  focal  length  in  meters. 

•  SubapertureDiameter  -  The  diameter  of  the  subaperture  in  meters. 

•  AverageWavelength  -  The  average  wavelength  in  meters. 

•  DetectorArraySize  -  The  square  root  of  the  number  of  pixels  in  each  subaperture  detector 
array.  The  CCD  array  is  assumed  square  with  no  dead  space  between  pixels. 

•  ReadNoise Variance  -  The  variance  for  the  additive  Gaussian  noise. 

•  AveragePhotonCount  -  The  average  photon  count  per  subaperture  for  the  sampling  period. 

•  WavefrontTilt  -  The  global  tilt  on  the  incident  wavefront  can  be  removed  to  simulate  the 
presence  of  a  tip-tilt  mirror  in  the  adaptive  optics  system. 

•  DetectorDutyCycle  -  The  duty  cycle  of  the  detector  array. 

•  SubapertureOverlap  -  The  intensity  can  spill  over  onto  adjacent  subapertures. 

The  vector  Xest  contains  the  estimated  centroid  shifts  for  all  subapertures  and  the  vector  Xact 
contains  the  actual  centroid  shifts  for  all  subapertures.  The  format  for  both  the  actual  centroid 
shift  vector  and  centroid  shift  estimate  vector  follows: 


S-octCl) 

msil) 

VactiX) 

my{l) 

®act  (2) 

m^(2) 

S/act(2) 

drnd  ^est  — 

mp{2) 

ms{M) 

3/oct(Af) 

mg{M) 
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where  a;act(l)  and  j/oct(l)  are  the  actual  x  and  y  shifts  in  meters  for  the  first  subaperture,  mg{l)  and 
mj7(l)  are  the  estimated  x  and  y  shifts  in  meters  for  the  first  subaperture,  M  is  the  total  number 
of  subapertures  in  the  wavefront  sensor  array,  and  the  subapertures  are  numbered  sequentially 
increasing  from  left  to  right,  top  to  bottom  as  depicted  in  Fig.  2.7. 
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CM 

CO 

D 

5 

CO 

7 

8 

9 

Figure  2.7  The  numbering  sequence  for  the  wavefront  sensor  array. 
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III.  Implementation  of  a  maximum  a-posteriori  estimator 

3.1  Introduction 

The  Hartmann  wavefront  sensor  estimates  the  centroid  shift  of  the  irradiance  distributions 
in  each  subaperture.  The  centroid  shift  information  is  then  used  to  estimate  the  wavefront  slope 
at  the  pupil  of  a  telescope  in  order  to  correct  for  wavefront  distortions  caused  by  atmospheric 
turbulence.  Centroid  shift  estimation  uses  a  moment  calculation  and  does  not  take  advantage  of 
the  correlation  properties  of  the  wavefront  slopes  over  the  subapertures  or  the  amount  of  light 
collected  by  the  wavefront  sensor.  A  MAP  estimator  derived  by  Scott  SaJlberg  [23]  incorporates 
statistical  knowledge  of  wavefront  slopes  and  light  levels  when  estimating  the  centroid  shifts.  The 
MAP  estimator  was  found  to  be  unbiased  and  the  MSE  performance  of  the  centroid  estimator  was 
its  upper  bound.  This  chapter  will  present  an  overview  of  the  MAP  estimator,  develop  the  equations 
necessary  to  implement  the  estimator  in  simulation,  and  evaluate  the  theoretical  performance  limits 
of  the  estimator. 

3.2  Maximum  a-posteriori  estimator 

The  MAP  estimator  derived  by  Sallberg  [23]  determines  the  MAP  shift  estimate  vector 
{xmap)  by  scaling  the  centroid  shift  estimate  vector  (xeat)  by  a  correction  factor  matrix  C~^: 

XmAP  ~  a  (3.1) 


where  the  correction  factor,  C  ^ ,  is  defined  [23]  as 


(3.2) 
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I  is  the  identity  matrix,  R  is  the  slope  correlation  matrix,  Cp  is  the  mean  square  spot  size,  and  K 


is  the  photon  count  matrix  defined  as 


Ki  0 


0  Kj 


(3.3) 


where  Kj  is  a  diagonal  matrix  containing  the  photon  count  in  the  jth  subaperture. 


S.3  Slope  correlation 

In  order  to  develop  a  simulation  for  the  MAP  estimator,  the  slope  correlation  matrix  R  must 
be  calculated.  The  goal  of  the  next  three  sections  is  to  derive  the  equations  needed  in  order  to 
calculate  slope  correlation  matrices  for  the  desired  atmospheric  statistics  and  subaperture  array 
size.  Referencing  Fig.  3.1,  a  general  solution  will  be  developed  to  find  slope  correlations  between 
any  two  subapertures  separated  by  a  distance  where  Xo  and  j/o  are  the  x  and  y 

distances  between  any  two  subapertures. 


Figure  3.1  The  slope  correlation  generalized  for  any  two  subapertures  separated  by  a  distance 
D.  The  coordinate  reference  is  centered  on  the  subapertures  with  positive  directions 
labeled. 
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The  slope  correlation  equations  are  developed  by  first  defining  the  slope  in  the  subaperture  as 
a  function  of  the  subaperture  weighting  function  Wi(x)  and  the  phase  in  the  subaperture  ^(x).  The 
wavefront  slope  for  the  ith  subaperture  can  be  shown  to  be  equivalent  to  computing  the  average 
wavefront  phase  gradient  over  the  wavefront  sensor  subaperture  [33]: 

Si  =  j  dS  Wi{S)  (V0(f)  •  di) ,  (3.4) 

where  di  is  a  unit  vector  designating  one  of  two  orthogonal  slope  components  for  the  ith  subaperture 
and  the  subaperture  weighting  function,  Wi{x),  is  normalized  such  that 

JdSWi{x)  =  l.  (3.5) 

Integrating  Eqn.  (3.4)  by  parts  [28]  allows  the  slope  to  be  expressed  in  terms  of  (j){x)  rather  than 
V</.(x): 

Si  =  -JdS  (VM^i(f)  •  di)0(f).  (3.6) 

The  slope  correlation  between  any  two  subapertures  i  and  j,  Rij,  is  the  expected  value  of  the 
product  of  the  slopes  in  the  two  subapertures  [22]: 

Rij  =  E\SiSj], 

=  ||df'dx(vW^i(f').dl)(viF,(f)-d;)E{.A(x')0(x)},  (3.7) 

where  s(i)  is  the  ith  subaperture  slope  measurement  in  the  di  direction  and  E[-]  is  the  expected 
value  operator. 
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Given  the  wavefront  phase  correlation  r^(5'  —  x)  =  E{(f)(x')(f>{x)}  and  the  definition  of  the 


structure  function  [11]: 

D^ix'  -  f)  =  2T^(x'  -x)-  2r^(0, 0), 


(3.8) 


the  slope  correlation  can  now  be  written  as 


Rij  =  J  j  d^'df  •  di)  (vw,(®)  •  4)  (r^(o,  0)  -  -  f)) . 


(3.9) 


Using  the  definition  of  the  subaperture  weighting  function  from  Eqn.  (2.4),  the  gradient  of 
the  weighting  function  was  found  to  be 

VM',(x,ri  =  i(i(x  +  |)-«(x-^))iect(2)x 

where  d  is  the  telescope  subaperture  diameter.  The  telescope  subaperture  diameter  is  used  in  this 
formulation  since  the  slope  correlation  values  will  be  calculated  across  the  aperture  of  the  telescope. 

Equations  to  calculate  the  slope  correlation  between  any  two  subapertures  can  now  be  devel¬ 
oped  using  Eqns.  (3.9)  and  (3.10).  It  should  be  noted  that  the  slope  correlation,  Rij,  is  dependent 
on  the  separation  of  the  ith  and  jth  subapertures,  Xg  and  yg,  and  the  slope  components,  di  and  dj, 
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in  each  subaperture.  Using  the  notation  Rij  =  R{xo,  Vo)  we  have, 


~  1 1 1 1 

(i  (*(*' + 1  -  - 1  -  •  '5) 

(*(*''  +  5  -  S'”)  -  ■'(*''  -  5  - 

X  (r,^(0. 0)  -  -x,y'  -  y)J . 


(3.11) 


S.4  Slope  correlation  matrix  development 

Before  specific  equations  can  be  developed  to  complete  the  correlation  matrix,  some  book¬ 
keeping  information  must  be  presented.  The  form  of  R  is  defined  by  the  ordering  of  the  elements 
of  8,  the  vector  containing  the  subaperture  wavefront  slopes: 


51 

52 

53 


S2M^-l 

52M2 


(3.12) 


where  si  is  the  x  directed  slope  for  the  1st  subaperture  (upper  left  in  Fig.  3.2)  and  S2  is  the 
y  directed  slope  for  the  1st  subaperture.  Since  the  array  of  subapertures  is  M  x  M,  with  each 
subaperture  having  an  x  and  y  directed  slope  measurement,  s  will  have  2M^  elements.  The  slope 
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correlation  matrix,  R  is  defined  as  The  slope  correlation  matrix  will  be  of  size  2M^  x  2M^ 

and  will  have  the  following  form: 


5iai 

S1S2 

S\S2M^ 

S2S1 

S2S2 

S2S2M^ 

S3S1 

S3S2 

SzS2M^ 

S4S1 

S4$2 

S432M2 

2M2-iai 

S2M^-1S2 

•  •  •  S2M^-IS2M^ 

S2M^Si 

S2M^S2 

S2M^S2M^ 

(3.13) 


Referring  to  Fig.  3.2  and  Eqn.  (3.13),  calculating  the  values  for  R  requires  the  development 
of  four  different  equations: 


1.  Correlation  between  two  subapertures  with  x  directed  slopes. 

2.  Correlation  between  two  subapertures  with  y  directed  slopes. 

3.  Correlation  between  subapertures  with  x  and  y  directed  slopes. 

4.  Correlation  between  subapertures  with  y  and  x  directed  slopes. 
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^^Subaperture 
*• — Full  Aperture 


Figure  3.2  Subaperture  slope  geometry.  The  numbering  depicts  orthogonal  components  of  the 
slope  measurements. 
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S.5  Slope  correlation  equations 


Equation  (3.11)  is  further  developed  for  specific  atmospheric  structure  functions  in  Appen¬ 
dices  A  and  B.  Appendix  A  develops  the  slope  correlation  using  Kolmogorov  statistics.  Appendix  B 
develops  the  slope  correlation  using  von  Karman  statistics.  The  following  is  a  summary  of  the  equa¬ 
tions  obtained  using  the  two  different  structure  functions. 


3.5.1  Kolmogorov  turbulence  statistics.  Using  the  Kolmogorov  structure  function  the 
correlation  between  x  directed  slopes  is 

RUxo,  %)  (^)  '  dA„(l  -  |A„  -  ilol) 

X  [2  {xl  +  A2)  "  -  ((f,  1)2  a2)  "  -  -  If  +  Al)  '] ,  (3.14) 

where  Xo  =  Xo/d  and  po  =  Vo/d  are  the  normalized  subaperture  separations  and  ro  is  the  Pried 

parameter.  The  correlation  between  y  directed  slopes  is 

Ryyi.Xo,%)  (^)  '  dA,(l  -  |A.  -  X,\) 

“  '<0'  Jxo-1 

X  [2 {yl  -H  A2)  '  -  ((j/„  -F  1)2  -h  A2)  "  -  ((y„  -  1)2  +  Al)  "] .  (3.15) 


The  correlation  between  x  and  y  directed  slopes  is 


RxyiXoi  Vo)  — ■ 


3.44/ d\t  p  fv^+t  ((  1\2  /  1  _  s2\6 

liL-i  '"'“*[((*'+2)  +(-2+*”-*) ) 

-  ((" + 5)' + (5 -’^)’) '-((»-  5)" + (-5 +*"-*)  7 
+((>'-l)'+(5+*‘'-*)T  ■ 


(3.16) 
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The  correlation  between  y  and  x  directed  slopes  is 


Ryxi^Ot  Vo) 


3.44/d\l  fi  ^1//  1^2  Z'  1  - 

liLi  j 

-((“'■*■  i)’  +  (I +»■>  - »)')  '-((*- 1)' + (-5 +»<■-»)’) 
+  ((“^-5)'  +  (5+»‘’-'')T]- 


(3.17) 


3.5.2  von  Karman  turbulence  statistics.  Similar  relationships  were  found  for  the  von 
Karman  structure  function.  The  correlation  between  x  directed  slopes  is 

R..(So, }.)  =0.08663(^)  ®  /’*”  -  |A,  -  fcl) 


rfv(^+i)^+(Aj,)^ 

Lc 

'  d,yiXo-iy+(Ay)^ 
Lo 


)*jtv.[2»^^5^3SS' 

)*jr,,.[2,ibZ(I^±<5E]], 


(3.18) 


where  Xg  =  Xg/d  and  yg  =  yg/d  are  the  normalized  subaperture  separations  and  Vg  is  the  Pried 
parameter.  The  correlation  between  y  directed  slopes; 


R„{ic,yo)  =0.08663(^)*^(EM)  r^”dA.(l-  |A.  -x„|) 

_  J5^  _ 


Lo 


=)!)V.,e[2xi^S±iS!]', 


(3.19) 
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The  correlation  between  x  and  y  directed  slopes  is 


R„{io, if.)  =0.08663(^)  ■  i (M)  I"'*’ 


+ 

+ 


-  V5/6 


■£'0 

I  L» 


The  correlation  between  y  and  x  directed  slopes  is 


R..(fe,fc)  =00S6e3(^)-l(M) 


2^d'v/(a!+|)=‘+(Hi/°-# 


'^-\)^+{-\+Vo-vy 


_  ^  _ 

_  ^<iA/(!c-^)^+(;+yo-i^^  ^ItTa/o 


(3.20) 


(3.21) 


3.5.3  Implementation.  The  correlation  between  subapertures  with  slopes  in  the  same 
direction  can  be  calculated  using  a  single  numerical  integration  where  as  the  correlation  between 
subapertures  with  x  and  y  directed  slopes  requires  two  dimensional  numerical  integration.  The 
correlation  matrix,  R,  for  the  Kolmogorov  structure  function,  once  generated,  can  simply  be  scaled 
as  the  Pried  parameter  (ro)  and  subaperture  diameter  (d)  changes.  However,  as  the  outer  scale  (Lg) 
and  subaperture  diameter  (d)  change  in  the  von  Karman  structure  function,  the  entire  correlation 
matrix  must  be  recalculated.  Also,  when  calculating  the  correlation  matrix  using  the  von  Karman 
structure  function,  a  series  expansion  to  the  second  order  of  the  modified  Bessel  function  of  the 
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second  kind  was  used  when  integrating  across  the  origin.  Without  a  series  expansion,  numerical 
integration  across  the  origin  results  in  an  infinite  value  which  is  a  result  of  being  equal  to 

infinity.  The  use  of  a  series  expansion  allowed  the  cancelation  of  the  term  that  caused  the  infinite 
value.  The  code  for  generating  the  correlation  matrices  is  listed  in  Appendix  E. 


3.6  Theoretical  MAP  estimator  performance 

In  order  to  verify  the  correlation  matrix  was  properly  developed,  the  theoretical  performance 
of  the  MAP  estimator  will  be  compared  to  the  values  obtained  by  Sallberg  [24].  The  behavior  of  the 
correction  factor  matrix  will  be  examined,  and  the  relative  mean  square  error  (MSE)  performance 
of  the  MAP  estimator  will  be  examined.  The  MSE  performance  will  be  evaluated  using  both 
Kolmogorov  statistics  and  von  Karman  statistics  in  the  correlation  matrix. 

3. 6. 1  Correction  factor  matrix.  First,  the  correction  factor  matrix,  will  be  examined. 
C~^  should  have  the  following  properties: 


0  as  if  0 
I  as  K  —*  oo. 


(3.22) 


As  the  average  photon  count  per  subaperture  approaches  zero,  the  data  cannot  be  trusted  and  the 
correction  factor  forces  the  MAP  estimate  to  zero.  As  the  average  photon  count  per  subaperture 
gets  very  large,  the  data  can  be  trusted  and  the  correction  factor  becomes  the  identity  matrix,  I, 
and  the  MAP  estimate  is  simply  the  centroid  estimate.  The  behavior  of  the  MAP  estimator  can 
be  verified  by  examining  the  behavior  of  the  correction  factor  matrix.  Mesh  plots  of  the  correction 
factor  matrix  as  defined  in  Eqn.  (3.2)  for  different  photon  counts  per  subaperture  {K)  are  shown 
in  Fig.  3.3.  As  the  average  photon  count  per  subaperture  increases,  the  correction  factor  becomes 
the  identity  matrix.  The  trend  verifies  the  expected  behavior  of  the  correction  factor  matrix. 
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Figure  3.3  Mesh  plots  of  the  correction  factor  matrix  for  various  values  of  K,  the  average  photon 

count  per  subaperture,  using  von  Karman  statistics  to  develop  the  correlation  matrix. 
The  subaperture  diameter  was  9.2cm,  the  Fried  coefficient  was  10cm,  and  the  outer 
scale  was  50m. 

3.6.2  MAP  estimator  performance  using  Kolmogorov  statistics.  Next,  the  performance 
MAP  estimator  will  be  examined  by  comparing  the  ratio  of  the  MSE  of  the  MAP  estimator  to 
the  MSE  of  the  centroid  estimator.  In  order  to  proceed,  relationships  for  the  MSE  of  the  centroid 
estimator  and  the  MAP  estimator  must  be  defined.  The  MSE  for  the  centroid  estimator  is  defined 
as  [24] 

MSE{x)  =  tT[ap^K-^],  (3.23) 
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and  the  MSE  for  the  MAP  estimator  is  defined  as  [24] 


MSE(f)  =  tr 


(3.24) 


where  tr[-]  is  the  trace  operator,  R  and  K  are  defined  by  Eqns.  (3.13)  and  (3.3),  and  (r|  is  the 
mean  square  spot  size. 


Figure  3.4  is  a  plot  of  the  error  ratio  as  a  function  of  average  subaperture  photon  count, 
K,  using  Kolmogorov  statistics  in  the  MAP  estimator  correlation  matrix.  The  correlation  matrix 


Figure  3.4  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture. 

Kolmogorov  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
d  =  9.2cm. 


was  computed  using  Kolmogorov  statistics  for  a  5  x  5  array  of  unobscured  subapertures  with 
subaperture  diameter  =  9.2cm.  The  various  curves  in  each  figure  correspond  to  a  particular  value 
of  the  ratio  of  the  mean  square  spot  motion,  to  mean  square  spot  size,  (Tp.  This  ratio  can  be 
expressed  as  a  ratio  of  subaperture  diameter,  d,  to  Pried  parameter.  To  [24], 


fi  _  ifmf  _ 


(3.25) 
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To  (meters) 

3 

0.0799 

2 

0.1018 

1 

0.1544 

1/2 

0.2340 

1/4 

0.3547 

Table  3.1  Pried  parameter  for  various  ratios  of  (T\la^  using  Kolmogorov  statistics  with  subper- 
ture  diameter  d  =  9.2cm. 

where  fi  is  the  subaperture  focal  length,  and 

=  12.803/(fcV®/3di/3)^  (3  26) 

is  the  measured  mean  square  angular  tilt  for  a  square  subaperture  [28].  The  factor  0.37  in  the 
denominator  of  Eqn.  (3.25)  is  a  result  of  matching  the  e~^  points  of  a  Gaussian  spot  distribution 
to  the  diffraction  pattern  of  a  square  subaperture  [18].  Table  3.1  lists  the  values  for  the  Pried 
parameter,  Tq,  associated  with  the  ratios  identifying  the  curves  in  Pig.  3.4. 

S.6.S  MAP  estimator  performance  using  von  Karman  statistics.  Pinally,  the  performance 
MAP  estimator  will  be  examined  by  comparing  the  ratio  of  the  MSE  of  the  MAP  estimator  to 
the  MSE  of  the  centroid  estimator  using  von  Karman  statistics  in  the  MAP  estimator  correlation 
matrix.  Pigure  3.5  is  a  plot  of  the  error  ratio  as  a  function  of  average  subaperture  photon  count,  K. 
The  correlation  matrix  was  computed  using  von  Karman  statistics  for  a  5  x  5  array  of  unobscured 
subapertures  with  outer  scale  Lg  —  50m  and  subaperture  diameter  d  =  9.2cm.  The  various  curves 
in  each  figure  correspond  to  a  particular  value  of  the  ratio  of  the  mean  square  spot  motion,  to 
mean  square  spot  size,  (Tp. 

A  relationship  must  now  be  found  to  relate  the  Pried  parameter,  rg,  to  when  von 

Karman  atmospheric  statistics  are  used.  The  first  step  is  to  write  Eqn.  (3.26)  as  a  function  of  the 
slope  correlation.  The  variable  JRxx  will  be  used  designate  the  normalized  mean  square  slope  and 
can  be  evaluated  for  Kolmogorov  statistics  using  Eqn.  (3.14)  when  the  normalized  offsets  {xg,yg) 
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Figure  3.5  Relative  MSB  performance  as  a  function  of  average  photon  count  per  subaperture. 

Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
outer  scale  Lg  =  50m  and  d  =  9.2cm. 


To  (meters) 

3 

0.0702 

2 

0.0896 

1 

0.1358 

1/2 

0.2058 

1/4 

0.3120 

Table  3.2  Pried  parameter  for  various  ratios  of  cr|/(rp  using  von  Karman  statistics  with  outer 
scale  Lo  =  50m  and  subperture  diameter  d  =  9.2cm. 

are  equal  to  zero,  To  =  1,  and  d  =  1.  Using  the  normalized  mean  square  slope,  Eqn.  (3.26)  can  be 
written  as. 


(3.27) 


where  2Rxx=12.803.  The  factor  of  2  is  required  since  the  slope  correlation  defined  by  Eqn.  (3.14) 
is  valid  for  slopes  in  a  single  orthogonal  direction  and  erg  is  defined  for  a  vector  parameter  in 
Eqn.  (3.26). 

A  similar  approach  can  be  used  to  derive  an  equation  for  cr|  using  von  Karman  statistics. 
The  variable  Rxx{d,Lo)  will  be  used  designate  the  normalized  mean  square  slope  for  a  specific  d 
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and  Lo  using  Eqn.  (3.18)  when  the  normalized  offsets  {xo,yo)  are  equal  to  zero.  The  mean  square 
angular  tilt  for  a  square  subaperture  in  Eqn.  (3.26)  can  be  written: 


2  _  ‘^Rxxjdi  Lq) 

‘  ~  evT 

The  factor  of  2  is  required  since  the  slope  correlation  defined  by  Eqn.  (3.18)  is  valid  for  slopes  in 
a  single  orthogonal  direction.  As  the  outer  scale  approaches  infinity,  the  value  of  the  normalized 
mean  square  slope  for  von  Karman  atmospheric  statistics  will  approach  the  value  of  the  normalized 
mean  square  slope  for  Kolmogorov  statistics: 

2R^x{d  =  l,Lo^oo)-*  12.803.  (3.29) 

The  ratio  between  mean  square  spot  motion  and  mean  square  spot  size  in  Eqn.  (3.25)  can  now 
be  written  in  terms  of  the  normalized  mean  square  slope  using  von  Karman  atmospheric  statistics: 


_  2Rxxid,  Lo)cP 
~  r®/^(0.37(27r))2' 


(3.30) 


For  a  given  ratio  of  the  Pried  parameter,  r^,,  can  be  determined.  Table  3.2  lists  the  values 

of  the  Pried  parameter.  To,  associated  with  the  ratios  identifying  the  curves  in  Fig.  3.5. 
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IV.  Verification  of  simulation  performance  and  test  plan 

4.1  Introduction 

This  chapter  verifies  the  accuracy  of  the  wavefront  sensor  simulation  and  maximum  a  posteri¬ 
ori  (MAP)  slope  estimator  implementation  as  well  as  outlining  a  test  plan  for  the  MAP  estimator. 
First,  the  wavefront  sensor  simulation  performance  was  examined.  The  accuracy  of  the  wavefront 
sensor  simulation  was  verified  by  comparing  the  shot  limited  performance  to  the  theoretical  bound. 
Since  the  wavefront  sensor  simulation  employs  a  centroid  estimator  to  determine  shift  estimates, 
the  performance  of  the  centroid  estimator  was  compared  to  the  shot  noise  limited  performance. 
Next,  the  effects  of  read  noise  and  atmospheric  turbulence  were  examined.  The  second  part  of 
this  chapter  verified  the  performance  of  the  MAP  estimator  implementation.  The  MAP  estimator 
mean  square  error  (MSE)  performance  was  evaluated  as  the  ratio  of  MAP  estimator  MSE  to  cen¬ 
troid  estimator  MSE.  The  curves  generated  from  the  simulations  were  compared  to  the  theoretical 
performance  curves.  The  final  part  of  this  chapter  outlines  a  test  plan  to  fully  evaluate  the  MAP 
estimator  performance  using  a  variety  of  wavefront  sensor  and  atmospheric  parameters. 

4.2  Hartmann  wavefront  sensor  simulation  verification 

In  order  to  verify  the  performance  of  the  wavefront  sensor  simulation,  the  shot  noise  limited 
performance  of  the  wavefront  sensor  was  examined.  The  shot  noise  limited  performance  was  evalu¬ 
ated  using  a  point  source  irradiance  while  ignoring  the  effects  of  read  noise  on  the  detection  process 
and  removing  any  effects  of  atmospheric  turbulence  from  the  incident  wavefront  phase.  Using  these 
assumptions,  each  subaperture  in  the  wavefront  sensor  array  forms  an  intensity  distribution  cen¬ 
tered  on  the  optical  axis.  While  the  intensity  distributions  are  centered  on  the  optical  axis  of  each 
subaperture,  shot  noise  in  the  detection  process  results  in  an  apparent  shift  off  the  optical  axis. 
The  shift  off  the  optical  axis  is  the  error  directly  resulting  from  the  presence  of  shot  noise  in  the 
detection  process.  This  error  is  the  shot  limited  performance  of  the  wavefront  sensor  simulation. 
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The  shot  noise  limited  performance  of  the  simulation  was  compared  to  the  theoretical  bounds 
in  order  to  verify  the  simulation  accuracy.  Since  the  shot  noise  limited  performance  of  a  2  x  2  CCD 
detector  array  (or  quad  cell  detector)  has  been  agreed  upon  by  several  sources  [17,27,31],  a  quad 
cell  detector  was  used  to  verify  the  simulation  performance.  The  theoretical  shot  noise  limited  root 
mean  square  error  (RMSE)  performance  for  a  subaperture  using  a  quad  cell  detector  in  the  absence 
of  atmospheric  turbulence  was  derived  by  Welsh  [31]  to  be 


where  K  is  the  average  photon  count  per  subaperture,  //  is  the  subaperture  focal  length,  A  is  the 
average  wavelength,  and  is  the  wavefront  sensor  subaperture  diameter.  The  RMSE  calculation 
in  Eqn.  (4.1)  is  the  Cramer-Rao  lower  bound  (CRLB)  for  any  unbiased  shift  estimator  [32].  In 
order  to  compare  the  shot  noise  limited  performance  of  the  simulation  to  the  theoretical  limit,  a 
centroid  shift  estimation  technique  developed  by  Tyler  [27]  was  used  in  place  of  the  first  moment 
calculation  defined  by  Eqn.  (2.12).  Tyler’s  method  optimizes  the  centroid  shift  estimation  for  a 
quad  cell  detector  by  scaling  the  ratio  of  the  difference  in  intensity  falling  on  either  half  of  the 
detector  (A/j,  and  AJy)  to  the  total  intensity  falling  of  the  detector  (/p)  by  a  constant  in  order  to 
obtain  the  centroid  shifts  in  the  x  and  y  directions.  The  centroid  shift,  x,  can  be  determined  for  a 
square  aperture  and  a  point  source  irradiance  using  the  relationship  developed  in  Appendix  C: 


A4M.  AlyfiX. 

T  j  ^  T  ^  y 


(meters), 


(4.2) 


where  //  is  the  subaperture  focal  length,  A  is  the  average  wavelength,  and  is  the  wavefront 
sensor  subaperture  diameter.  It  should  be  noted  that  this  shift  estimator  was  derived  for  a  quad 
cell  detector  with  shot  noise  and  is  not  valid  when  read  noise  and  atmospheric  turbulence  are 
included  in  the  wavefront  sensor  simulation. 
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Equation  (4.2)  was  used  to  calculate  the  centroid  shift  in  the  wavefront  sensor  simulation.  The 
simulation  RMSE  performance  was  compared  to  the  theoretical  RMSE  performance.  Figure  4.1 
shows  the  wavefront  sensor  simulation  performance  matched  the  shot  noise  limited  performance. 


The  results  were  obtained  by  averaging  the  RMSE  for  500  simulation  runs  at  each  value  of  K. 


Figure  4.1  Actual  vs  Theoretical  RMSE  for  a  wavefront  sensor  using  a  quad  cell  detector.  The 
wavefront  sensor  simulation  used  a  modified  shift  estimator  and  was  run  in  the  absence 
of  read  noise  or  atmospheric  turbulence. 


4-3  Centroid  estimator  RMSE  performance 

4-3.1  Shot  noise  effects.  The  performance  of  the  centroid  estimator  was  examined  next. 
As  seen  in  Fig.  4.1,  the  RMSE  performance  of  the  wavefront  sensor  with  quad  cell  detector  ap¬ 
proached  the  Cramer-Rao  bound  when  a  modified  shift  estimator  was  used  to  find  the  centroid 
shift.  The  wavefront  sensor  simulation  employed  a  centroid  estimator  that  uses  a  moment  calcu¬ 
lation  described  by  Eqn.  (2.12)  to  find  the  centroid  shifts.  Since  the  centroid  estimator  was  not 
optimized  for  shot  noise  limited  performance,  as  was  the  estimator  derived  in  Appendix  C,  the 
use  of  a  centroid  estimator  should  result  in  shift  estimates  with  more  error.  At  best,  the  centroid 
estimator  performance  should  match  shot  noise  limited  performance  as  the  number  of  pixels  in  the 
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detector  array  gets  very  large.  Figure  4.2  shows  the  performance  of  the  centroid  estimator  as  well 
as  the  shot  noise  limited  performance  developed  in  the  previous  section. 


Figure  4.2  Comparison  of  the  centroid  estimator  and  the  Cramer-Rao  RMSE  performance  for 
a  wavefront  sensor  using  a  quad  cell  detector.  The  wavefront  sensor  was  a  5  x  5 
unobscured  array  of  detectors  without  read  noise  or  atmospheric  turbulence. 

4-3.2  Read  noise  effects.  Since  read  noise  is  present  in  the  detection  process,  the  effects 
it  has  on  both  the  theoretical  bounds  and  centroid  estimator  performance  were  examined.  The 
addition  of  read  noise  in  the  detection  process  should  increase  the  shift  estimate  error.  First,  the 
effects  on  the  theoretical  bounds  will  be  examined.  If  the  total  error  in  the  detection  process  is 
assumed  to  result  from  shot  and  read  noise,  the  MSE  can  be  written  as  a  sum  of  the  MSE  due  to 
shot  noise  and  the  MSE  due  to  read  noise  [17]: 

+  <^rn  (4.3) 

Parenti  developed  a  pair  of  equations  to  describe  the  error  resulting  from  shot  and  read  noise  when 
the  wavefront  sensor  used  a  quad  cell  detector  in  each  subaperture  [17].  The  MSE  resulting  from 
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shot  noise  is, 
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and  the  MSE  resulting  from  read  noise  is, 


OO 

r  ^ 

K2  > 

\2'Kdw} 

(meters^). 


(4.5) 


where  K  is  the  average  photon  count,  is  the  variance  of  the  read  noise  for  a  single  detector 
pixel,  fi  is  the  subaperture  focal  length,  is  the  wavefront  sensor  subaperture  diameter,  A  is  the 
average  wavelength,  and  To  is  the  Pried  parameter.  If  no  turbulence  is  assumed  (ro  — >  oo),  the 
expression  for  the  MSE  in  Eqn.  (4.4)  is  consistent  with  Eqn.  (4.1). 

Using  Eqns.  (4.4)  and  (4.5),  the  theoretical  performance  of  a  wavefront  sensor  using  a  quad 
cell  detector  in  each  subaperture  was  examined  under  two  conditions.  First,  the  performance  was 
examined  as  read  noise  increased  in  the  absence  of  atmospheric  turbulence  (fp  — >  oo).  Then,  the 
performance  was  examined  as  the  level  of  atmospheric  turbulence  increased  in  the  absence  of  read 
noise  {(tr  =  0). 

The  curves  in  Fig.  4.3  show  the  effects  of  increasing  read  noise  variance  on  theoretical  RMSE 
performance.  At  low  photon  counts,  the  read  noise  significantly  increases  the  shift  error.  As 
the  photon  count  increases,  the  MSE  due  to  read  noise  decreases  by  a  factor  of  K^,  and  the 
overall  MSE  performance  converges  to  the  shot  noise  limited  performance  of  the  centroid  estimator. 
Similar  results  were  obtained  with  the  wavefront  sensor  simulation.  The  performance  of  the  centroid 
estimator  for  increasing  levels  of  read  noise  variance  is  shown  in  Fig.  4.4.  At  low  photon  counts,  the 
read  noise  significantly  effects  the  performance  of  the  centroid  estimator  and  as  the  photon  count 
increases,  the  error  curves  converge  to  the  shot  noise  limited  performance  curve  of  the  centroid 
estimator  (ct^  =  0). 


4-5 


0  20  40  60  80  100  120  140  160  180  200 


Average  photon  count  per  subaperture,  K 

Figure  4.3  Theoretical  RMSE  performance  as  read  noise  variance  {(t%)  increases  for  a  quad  cell 
detector.  The  wavefront  sensor  was  a  5  x  5  unobscured  array  of  subapertures  without 
atmospheric  turbulence. 

4-3.3  Atmospheric  turbulence  effects.  The  theoretical  and  centroid  estimator  RMSE  per¬ 
formance  in  the  presence  of  atmospheric  turbulence  were  examined  next.  Atmospheric  turbulence 
deforms  the  overall  shape  of  the  intensity  distribution  that  forms  on  the  detector  array.  In  addition 
to  the  shot  noise,  the  deformation  contributes  to  the  shift  error  during  the  detection  process.  The 
curves  in  Fig.  4.5  show  the  effects  of  increasing  atmospheric  turbulence  on  theoretical  RMSE  per¬ 
formance.  The  individual  curves  correspond  to  a  particular  ratio  of  crl/cTp.  The  theoretical  RMSE 
curves  were  calculated  using  Eqn.  (4.4)  using  values  of  r©  found  in  Table  3.2.  The  theoretical  RMSE 
increased  as  the  atmospheric  turbulence  increased.  Similar  results  were  obtained  with  the  wave- 
front  sensor  simulation.  The  curves  in  Fig.  4.6  show  the  centroid  estimator  RMSE  performance  as 
atmospheric  turbulence  increases. 

4-4  MAP  estimator  performance 

4.4-i  Introduction.  In  order  to  verify  the  MAP  estimator  performance,  the  combination  of 
the  phase  screen  generation  program,  the  wavefront  sensor,  and  the  MAP  estimator  implementation 
was  used  to  generate  a  set  of  baseline  performance  curves.  The  baseline  performance  curves  were 
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Figure  4.4  RMSE  performance  for  a  centroid  estimator  using  a  quad  cell  detector  with  read  noise. 

The  curves  correspond  to  a  particular  read  noise  variance  (cfl).  As  increases,  the 
RMSE  increases.  The  wavefront  sensor  was  a  5  x  5  unobscured  array  of  subapertures 
without  atmospheric  turbulence. 

obtained  using  a  set  of  parameters  that  matched  the  theoretical  assumptions  used  in  the  derivation 
of  the  MAP  estimator.  The  curves  generated  in  simulation  should  match  the  theoretical  curves 
shown  in  Fig  3.5. 

4-4'2  Simulation  flowchart.  Figure  4.9  is  a  flow  chart  for  the  shell  program  that  imple¬ 
ments  the  MAP  estimator  in  simulation.  After  choosing  a  desired  ratio  of  «rg/cr|,  rg  was  determined 
using  Eqn.  (3.30).  Then,  a  series  of  phase  screens  were  generated,  and  the  correlation  matrix  was 
calculated  for  the  MAP  estimator  implementation.  For  each  phase  screen,  a  centroid  shift  estimate 
vector  (xesi)  and  an  actual  centroid  shift  vector  (»oc<)  was  determined  using  the  wavefront  sensor 
simulation.  The  centroid  shift  estimate  vector  was  then  used  to  calculate  the  MAP  shift  estimate 
vector  (xmap)-  The  average  MSE  for  the  centroid  and  MAP  estimates  were  then  calculated  using 
(xact)  as  the  correct  shift  vector.  This  process  was  repeated  for  values  of  K  ranging  from  1  to  200 
photons  per  subaperture. 
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Figure  4.5  RMSE  performance  for  a  quad  cell  detector  array  with  no  read  noise.  The  curves  cor¬ 
respond  to  a  particular  value  of  tr^/tTp.  The  wavefront  sensor  was  a  5  x  5  unobscured 
array  of  subapertures. 

4-4-3  Simulation  performance.  In  theory,  the  ratio  of  the  MAP  estimator  MSE  to  the 
centroid  estimator  MSE  should  match  the  theoretical  curves  shown  in  Fig.  3.5.  The  parameters 
used  to  generate  the  baseline  performance  curves  are  listed  in  Table  4.1.  The  physical  wavefront 
sensor  parameters  were  chosen  to  match  the  Generation  III  Hartmann  wavefront  sensor  parameters 
used  at  the  Starfire  Optical  Range  [19].  However,  instead  of  using  a  16  x  16  array  of  subapertures, 
the  simulation  used  a  5  x  5  array  in  order  to  keep  the  array  sizes  manageable.  (Larger  array  sizes 
would  have  exceeded  the  resident  memory  capabilities  of  the  computational  platforms  used,  greatly 
increasing  processing  time  for  each  simulation.)  In  addition  to  the  decreased  sensor  array  size, 
certain  parameters  were  changed  to  approximate  the  theoretical  assumptions  made  by  Sallberg  [23]: 
Infinite  resolution  of  the  subaperture  detector  array  was  approximated  using  a  20  x  20  detector 
array,  the  read  noise  was  eliminated,  and  overlap  between  subapertures  was  not  allowed. 

The  baseline  performance  curves  of  the  ratio  of  the  MAP  estimator  MSE  to  the  the  centroid 
estimator  MSE  as  a  function  of  average  photon  count  per  subaperture,  K  are  shown  in  Fig.  4.7.  The 
dashed  curves  correspond  to  the  theoretical  performance  shown  in  Fig.  3.5  while  the  solid  curves 
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Figure  4.6  RMSE  performance  for  a  centroid  estimator  using  a  quad  cell  detector  array  with 
no  read  noise.  The  curves  correspond  to  a  particular  value  of  tr^/o-p.  The  wavefront 
sensor  was  a  5  x  5  unobscured  array  of  subapertures. 

correspond  to  the  simulation  performance.  Von  Karman  statistics  were  used  for  a  5  x  5  array  of 
subapertures  with  Lo  =  50m  and  d  =  9.2cm.  The  pairs  of  curves  correspond  to  a  particular  ratio 
oiallal. 

The  MAP  estimator  MSE  performance  was  upper  bounded  by  the  centroid  estimator  MSE 
as  predicted  by  Sallberg  [23].  While  the  baseline  performance  curves  generated  by  the  simulation 
should  match  the  theoretical  curves,  the  expected  relative  MAP  estimator  MSE  performance  was 
not  quite  achieved.  In  fact,  the  relative  MAP  estimator  MSE  performance  in  simulation  appeared 
to  perform  better  than  the  theoretical  expectations.  There  were  several  reasons  why  the  curves 
generated  from  simulation  did  not  fall  directly  on  top  of  the  theoretical  curves. 

The  first  reeison  the  theoretical  and  simulation  curves  do  not  exactly  match  was  that  the 
actual  value  of  mean  square  motion,  cr^,  to  the  mean  square  spot  size,  cr^,  calculated  from  the 
phase  screens  did  not  match  the  theoretical  ratio.  The  value  for  the  Fried  parameter,  r*o,  was 
determined  from  the  ratio  of  rising  Eqn.  (3.30).  The  same  To  was  used  both  in  the  phase 
screen  generation  and  in  the  calculation  of  the  correlation  matrix  for  the  MAP  estimator.  Using 
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Phase  screen  generator 

Subaperture  diameter  (d) 
Outer  scale  (Lq) 

Fried  Parameter  (to) 

9.2cm 

50m 

determined  by 

Wavefront  sensor 

Subaperture  array  size 

5x5 

Subaperture  width  (dw) 

200E-6m 

Subaperture  focal  length  (/j) 

650E-5m 

Average  wavelength  (A) 

550nm 

Subaperture  detector  array 

20  X  20 

CCD  duty  cycle 

100  percent 

Read  noise  (cr^) 

Variance  =  0 

Wavefront  tilt 

Included 

Subaperture  overlap 

Not  allowed 

MAP  estimator 

Atmospheric  statistics 

von  Karman 

Outer  scale  (£<>) 

50m 

Pried  parameter  (to) 

determined  by  cr^  /cr^ 

Table  4.1  Parameters  used  to  verify  the  simulations  match  the  theoretical  performance  limits. 


the  phase  screens  generated  from  the  value  of  To,  the  actual  mean  square  spot  motion  and  mean 
square  spot  size  were  measured  and  the  actual  ratio  of  (rl/a^  was  calculated.  The  theoretical  ratios 
did  not  always  match  the  calculated  ratio. 

The  second  reason  the  theoretical  and  simulation  curves  do  not  exactly  match  was  that 
the  MAP  estimator  was  not  necessarily  a  minimum  mean  square  error  (MMSE)  estimator.  A 
MAP  estimator  that  uses  parameters  with  jointly  Gaussian  distributions  will  result  in  a  MMSE 
estimator  [15].  However,  the  development  of  the  MAP  estimator  did  not  use  jointly  Gaussian 
distributions.  Shot  noise  was  modeled  with  a  Poisson  process  and  the  irradiance  distribution  was 
only  approximated  by  a  Gaussian  distribution.  Since  the  development  of  the  MAP  estimator  did 
not  use  jointly  Gaussian  distributions,  the  theoretical  performance  defined  by  Eqn.  (3.24)  was  not 
necessarily  the  MMSE  performance  for  the  MAP  estimator.  As  a  result,  the  actual  MSE  for  the 
MAP  estimator  in  simulation  may  fall  above  or  below  the  theoretical  MSE  performance  curve. 

Finally,  there  were  errors  associated  with  the  discrete  nature  of  the  simulations.  The  sampling 
of  the  phase  screen  was  done  at  the  minimum  Nyquist  rate  [30].  The  wavefront  sensor  used  a  discrete 
Fourier  Transform  to  find  the  intensity  distributions  for  each  subaperture.  The  calculation  of  the 
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Figure  4.7  Ratio  of  MSB’s  for  the  MAP  and  the  centroid  estimator  as  a  function  of  average 
photon  count  per  subaperture,  K.  Von  Karman  statistics  were  used  for  a  5  x  5  array 
of  subapertures  with  Lo  =  50m  and  d  =  9.2cm.  The  curves  correspond  to  a  particular 
ratio  of 

correlation  values  for  von  Karman  atmospheric  statistics  included  a  finite  series  approximation  of 
a  modified  Bessel  function  of  the  second  kind  as  well  as  numerical  integration.  While  steps  were 
taken  to  minimize  the  errors  associated  with  each  of  these  simulations,  it  must  be  recognized  that 
errors  exist  and  propagate  through  the  simulations. 

4-5  Test  plan 

Two  formats  were  used  to  present  the  simulation  results.  Absolute  error  plots  will  show 
the  MAP  and  centroid  estimator  MSB  as  a  function  of  average  photon  count  per  subaperture,  K. 
Relative  error  plots  will  show  the  ratio  of  the  MAP  estimator  MSB  to  the  centroid  estimator  MSB 
as  a  function  of  K.  An  improvement  in  the  relative  MAP  estimator  MSB  performance  does  not 
indicate  an  absolute  improvement  in  the  MAP  estimator  MSB  performance.  It  simply  means  the 
the  MAP  estimator  MSB  performance  relative  to  the  centroid  estimator  MSB  performance  has 
improved.  It  is  possible  for  the  relative  MAP  estimator  MSB  performance  to  improve  while  the 
absolute  MAP  estimator  MSB  performance  actually  gets  worse. 
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The  goal  of  the  remaining  simulations  was  to  examine  the  effects  of  changing  the  wavefront 
sensor  parameters  listed  in  Table  4.1  on  the  MSE  performance.  The  eventual  goal  was  to  test  the 
MSE  performance  using  actual  wavefront  sensor  and  atmospheric  parameters.  The  following  list 
shows  the  order  in  which  the  wavefront  sensor  and  atmospheric  parameters  will  be  examined  in 
simulation. 

•  Remove  global  tilt  on  the  incident  wavefront.  This  was  done  to  simulate  the  presence  of  a 
tip-tilt  mirror  in  the  adaptive  optics  system. 

•  Allow  overlap.  The  effects  of  allowing  light  to  spill  over  onto  adjacent  subapertures  were 
examined.  In  actual  wavefront  sensors,  there  is  no  way  to  prevent  spill  over  from  occurring 
on  the  detector  array  short  of  turning  off  the  pixels  bordering  adjacent  subapertures. 

•  Decrease  the  size  of  the  detector  array.  The  MAP  estimator  development  assumed  the  sub¬ 
aperture  detector  array  had  infinite  resolution.  The  infinite  array  was  approximated  using  a 
20  X  20  pixel  array  in  each  subaperture.  However,  typical  wavefront  sensors  utilize  quad  cell 
detectors  or  4  x  4  detector  arrays  [9,19].  The  simulations  will  investigate  the  performance 
for  8  X  8,  4  X  4,  and  2x2  detector  arrays  in  each  subaperture. 

•  Add  read  noise.  Read  noise  is  a  source  of  error  for  all  CCD  detection  devices.  Using  the 
4x4  and  2x2  detector  arrays,  the  effects  of  read  noise  was  examined  for  variances  of  2,  5, 
and  10  (photons/pixel)^. 

•  Incorrect  atmospheric  parameters  in  the  MAP  estimator.  Up  until  this  point  in  time,  it 
was  assumed  that  the  atmospheric  parameters  where  known  with  a  high  degree  of  certainty. 
While  the  wavefront  sensor  can  be  used  to  estimate  atmospheric  parameters  [3,6,25],  error 
within  the  wavefront  sensor  limits  the  degree  of  certainty  to  which  the  data  can  be  trusted. 
Therefore,  the  MAP  estimator  was  implemented  using  incorrect  atmospheric  parameters. 

•  Modified  detector  array.  A  modified  4x4  detector  array  was  tested.  The  detector  was 
designed  such  that  the  four  center  pixels  were  active  and  the  pixels  around  the  border  of  the 
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array  were  turned  off.  The  goal  with  this  design  was  to  minimize  electronic  crosstalk  between 
pixels  in  adjacent  subapertures.  This  configuration  was  used  in  the  SOR  Generation  III 
tests  [19]. 


Figure  4.8  The  figure  on  the  right  is  the  standard  4x4  detector  array  while  the  figure  on  the 
left  is  the  modified  SOR  Generation  III  detector. 
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Figure  4.9  A  flow  chart  for  the  combined  operation  of  the  phase  screen  generation  program,  the 
wavefront  sensor  simulation,  and  the  MAP  estimator  implementation. 
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V.  Results  and  Analysis 


5.1  Introduction 

To  this  point,  the  performance  of  the  MAP  slope  estimator  has  been  validated  using  wavefront 
sensor  parameters  that  approximated  the  assumptions  made  in  the  derivation  of  the  MAP  estimator: 
The  detection  array  was  assumed  large,  read  noise  was  not  included  in  the  detection  of  intensity 
distributions,  overlap  between  subapertures  was  prohibited,  and  the  atmospheric  parameters  were 
known  with  a  high  degree  of  certainty.  This  chapter  examined  the  MAP  and  centroid  estimator 
MSE  performance  as  the  idealized  wavefront  sensor  parameters  were  changed  to  reflect  parameters 
that  could  be  found  in  realizable  wavefront  sensors.  First,  the  global  tilt  was  removed  from  the 
incident  wavefront  to  simulate  the  presence  of  a  tip-tilt  mirror  in  the  adaptive  optics  system.  Next, 
the  intensity  distributions  were  allowed  to  overlap  into  adjacent  subapertures.  The  number  of  pixels 
in  each  subaperture  detector  array  was  decreased.  Finally,  read  noise  was  added  to  the  detection 
process.  In  addition  to  changing  the  wavefront  sensor  parameters,  the  atmospheric  parameters 
were  changed  as  well.  Up  to  this  point,  it  was  assumed  that  the  atmospheric  parameters  were 
known  with  a  high  degree  of  certainty.  Since  this  is  rarely  true,  the  effects  of  using  incorrect 
Fried  parameters  and  outer  scales  in  the  MAP  estimator  implementation  were  examined.  The  last 
series  of  simulations  examined  the  performance  of  the  MAP  estimator  when  a  modified  subaperture 
detector  array  was  used. 

5.2  Global  tilt 

Since  adaptive  optics  systems  typically  employ  a  tip-tilt  mirror  [13],  the  effect  of  removing 
the  global  tilt  from  the  incident  wavefront  was  examined.  The  theoretical  results  presented  by 
Sallberg  [24]  showed  that  the  relative  MAP  estimator  MSE  performance  decreased  when  the  global 
tilt  was  removed  from  the  incident  wavefront.  The  decreased  relative  MAP  estimator  MSE  per¬ 
formance  resulted  from  a  loss  of  slope  correlation  between  subapertures  due  to  the  tilt  removal. 
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Similar  results  should  occur  in  simulation.  Removing  the  global  tilt  from  the  incident  wavefronts 
should  result  in  a  decrease  in  the  relative  MAP  estimator  MSE  performance. 


Figure  5.1  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  for  tilt  removed  wavefronts.  The  curves  correspond  to  a  particular  ratio  of 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  subapertures  with  Lq  =  50m 
and  d  =  9.2cm. 

Figure  5.1  shows  the  relative  MSE  performance  when  the  parameters  listed  in  Table  4.1  were 
used  with  global  tilt  removed  from  the  incident  wavefronts.  There  was  no  significant  decrease  in 
relative  MAP  estimator  MSE  performance  when  the  global  tilt  was  removed  from  the  incident 
wavefronts.  The  results  presented  by  Sallberg  [24]  used  Kolmogorov  atmospheric  statistics.  Since 
von  Karman  atmospheric  statistics  were  used  in  the  simulations,  there  was  less  slope  correlation 
between  subapertures.  The  decrease  in  slope  correlation  resulting  from  the  removal  of  global  tilt 
from  the  incident  wavefronts  had  minimal  effect  on  the  relative  MAP  estimator  MSE  performance. 

5.3  Intensity  overlap 

The  effect  of  allowing  the  intensity  distributions  to  spill  over  onto  adjacent  subaperture  de¬ 
tector  arrays  was  examined  next.  While  the  ratio  of  mean  square  spot  motion  to  mean  square  spot 
size  has  been  used  to  identify  specific  performance  curves,  the  actual  wavefront  sensor  parameters 
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will  determine  the  actual  mean  square  spot  size  and  mean  square  spot  motion.  If  a  wavefront  sensor 
is  correctly  designed  to  operate  below  a  known  maximum  level  of  turbulence,  the  mean  square  spot 
motion  within  a  particular  subaperture  will  be  small  enough  to  prevent  the  intensity  distribution 
from  overlapping  onto  adjacent  subaperture  detector  arrays.  Since  the  wavefront  sensor  parame¬ 
ters  used  for  these  simulations  were  optimized  for  the  selected  turbulence  levels,  there  should  be 
no  significant  decrease  in  relative  MAP  estimator  MSE  performance  due  to  overlap. 


Figure  5.2  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture 
count,  K,  for  overlap  allowed.  The  curves  correspond  to  a  particular  value  of  the 
ratio  of  Von  Karman  statistics  were  used  for  a  5  x  5  array  of  subapertures 

with  Lo  =  50m  and  d  =  9.2cm. 

Figure  5.2  shows  the  relative  MAP  estimator  MSE  performance  when  the  parameters  listed 
in  Table  4.1  were  used  while  allowing  intensity  distributions  to  overlap  onto  adjacent  subapertures. 
There  was  only  a  slight  decrease  in  relative  MAP  estimator  MSE  performance.  The  SOR  Generation 
III  parameters  resulted  in  a  spot  that  was  about  1/6  the  width  of  the  subaperture  diameter.  The 
largest  value  of  cr^/cTp  used  in  the  simulations  constrained  the  mean  square  spot  motion  to  the 
center  25  percent  of  the  detector  array  99  percent  of  the  time.  As  a  result  of  the  limited  mean 
square  spot  motion  on  the  detector  array,  overlap  had  little  effect  of  the  relative  MAP  estimator 
MSE  performance.  While  the  relative  MAP  estimator  MSE  performance  showed  minimal  change. 
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there  was  a  slight  increase  in  the  absolute  MSE  performance  of  the  MAP  estimator  when  overlap 
vas  allowed.  Figure  5.3  shows  the  increase  in  MAP  estimator  MSE  performance  when  overlap  was 
allowed  for  cr^lo-p  =  1. 


Figure  5.3  MAP  estimator  MSE  performance  as  a  function  of  average  photon  count  per  subaper¬ 
ture  count,  K,  for  o’c/^'p  =  1-  Karman  statistics  were  used  for  a  5  x  5  array  of 
unobscured  subapertures  with  Lq  =  50m  and  d  =  9.2cm. 

In  order  to  observe  a  decrease  in  the  relative  MAP  estimator  MSE  performance,  the  wavefront 
sensor  parameters  were  changed  to  be  non-optimal  for  the  selected  turbulence  levels  by  increasing 
the  subaperture  focal  length  by  a  factor  of  five.  This  will  cause  the  irradiance  distribution  to 
overlap  onto  adjacent  subapertures  as  it  roams  about  the  detector  array.  For  the  same  ratio  of 
<Tg/(Tp  =  1  used  in  the  previous  simulation,  both  the  mean  square  spot  size  and  the  mean  square 
spot  motion  were  larger  due  to  the  increased  subaperture  focal  length.  The  effect  of  overlap  on  the 
relative  MAP  estimator  MSE  performance  became  obvious.  Figure  5.4  shows  that  when  the  focal 
length  was  increased  to  reduce  the  optimality  of  the  wavefront  sensor  design,  the  intensity  overlap 
into  adjacent  subapertures  significantly  decreased  the  relative  MAP  estimator  MSE  performance 
as  expected. 
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Figure  5.4  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture 
count,  K.  The  subaperture  focal  length  has  been  increased  to  150  x  subaperture 
diameter.  Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subaper¬ 
tures  with  (Tc/(Tp  =  1,  Lg  =  50m  and  d  =  9.2cm. 

5.4  Detector  array  size 

For  the  derivation  of  the  MAP  estimator,  the  detection  array  in  each  subaperture  was  assumed 
to  have  an  infinite  resolution.  Wavefront  sensors  typically  use  a  small  number  of  CCD  detector  pixels 
in  each  subaperture.  The  number  of  CCD  pixels  per  subaperture  is  limited  by  many  physical  and 
wavefront  sensor  design  related  factors  including  the  physical  size  of  each  CCD  pixel,  how  closely 
the  CCD  pixels  can  be  packed  together,  and  the  read  noise  variance  for  each  pixel.  Wavefront 
sensors  typically  use  a  quadcell  or  4  x  4  subaperture  detector  array. 

The  centroid  estimator  MSE  performance  should  increase  as  the  number  of  pixels  in  each 
subaperture  detector  array  decreases.  As  the  error  in  the  centroid  shift  estimates  increases  the 
MAP  shift  estimates  will  not  be  as  accurate.  The  relative  MAP  estimator  MSE  performance 
should  decrease  as  the  number  of  pixels  in  each  subaperture  detection  array  decreases. 

5.4.1  Relative  MSE  for  a  fixed  ratio  of  a^/up.  The  next  series  of  simulations  examined 
the  relative  MAP  estimator  MSE  performance  using  parameters  listed  in  Table  4.1  with  global 
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tilt  removed,  overlap  allowed,  and  subaperture  detection  array  sizes  of  8  x  8,  4  x  4,  and  2x2. 
Figures  5.5  -  5.8  show  the  relative  MAP  estimator  MSE  performance  as  a  function  of  average 
photon  count  per  subaperture,  K,  for  a  specific  value  of  a^lcTp,  with  tilt  removed  wavefronts  and 
overlap  allowed.  The  individual  curves  in  each  figure  correspond  to  the  labeled  subaperture  detector 
array  size.  For  all  levels  of  atmospheric  turbulence,  the  relative  MAP  estimator  MSE  performance 
decreased  as  the  detector  array  size  decreased.  The  grouping  of  the  20  x  20  and  8x8  curves 
and  the  4x4  and  2x2  curves  resulted  from  increased  centroid  shift  estimate  error.  The  denser 
20  X  20  and  8x8  detector  arrays  allowed  for  a  significantly  lower  centroid  shift  estimate  error 
than  the  4x4  and  2x2  detector  arrays.  The  increased  shift  error  decreased  the  relative  MSE 
performance  for  the  MAP  estimator. 


Figure  5.5  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  for  cr\l(Tp  =  2.  The  curves  correspond  to  a  specific  CCD  detector  array  size.  Von 
Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with  Lo  = 
50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 


5.4-2  Absolute  MSE  for  a  fixed  ratio  ofa^/ap.  The  absolute  MAP  and  centroid  estimator 
MSE  performance  using  the  ratio  of  crj/aj  =  1/2  with  tilt  removed  wavefronts  and  overlap  allowed 
is  shown  in  Fig.  5.9.  As  the  number  of  pixels  in  the  detector  array  decreased,  both  the  MAP  and 
centroid  estimator  MSE  performance  increased  as  expected.  While  the  relative  MAP  estimator 
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Figure  5.6  Relative  MSB  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  for  (r'^ldp  =  1.  The  curves  correspond  to  a  specific  CCD  detector  array  size. 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 


MSB  performance  in  Fig.  5.7  indicated  a  slight  decrease  as  the  number  of  pixels  in  the  detector 
array  decreased  for  =  1/2,  the  absolute  MAP  and  centroid  estimator  MSB  performance 

increased  by  an  order  of  10^  between  the  2x2  and  the  20  x  20  CCD  array.  This  shows  that 
the  resolution  of  the  detector  array  contributes  significantly  to  the  accuracy  of  both  the  MAP  and 
centroid  shift  estimates. 


5.4-3  Relative  MSB  for  a  fixed  detector  array  size.  The  same  MSB  performance  data  was 
presented  in  another  format.  Figures  5.10  -  5.12  show  the  relative  MAP  estimator  MSB  performance 
as  a  function  of  average  photon  count  per  subaperture,  K,  for  a  specific  detector  array  size  with  tilt 
removed  wavefronts  and  overlap  allowed.  The  curves  correspond  to  a  specific  value  of  (r^/cTp.  As 
expected,  the  relative  MAP  shift  estimator  MSB  performance  decreased  as  the  level  of  atmospheric 
turbulence  increased.  For  detector  arrays  with  few  pixels,  there  was  only  a  slight  decrease  in  relative 
MAP  estimator  MSB  performance  as  the  turbulence  levels  increased. 
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Average  photon  count  per  subaperture,  K 


Figure  5.7  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  for  al/cTp  =  1/2.  The  curves  correspond  to  a  specific  CCD  detector  array  size. 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 


5.4-4  Absolute  MSE  for  a  fixed  detector  array  size.  The  absolute  MSE  performance  for  the 
MAP  and  centroid  estimators  using  a  4  x  4  detector  array  with  tilt  removed  wavefronts  and  overlap 
allowed  for  various  turbulence  levels  is  shown  in  Fig.  5.13.  As  the  level  of  atmospheric  turbulence 
increased,  the  MAP  and  centroid  estimator  MSE  performance  increased  as  expected.  While  the 
relative  MAP  estimator  MSE  performance  only  showed  a  slight  change  for  differing  atmospheric 
turbulence  levels  in  Fig.  5.11,  the  absolute  MAP  and  centroid  estimator  MSE  performance  increased 
by  an  order  of  10  as  the  ratio  of  cr^/cTp  increased  from  1/4  to  2. 
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Figure  5.8  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  for  (Tg/(Tp  =  1/4.  The  curves  correspond  to  a  specific  CCD  detector  array  size. 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 


Figure  5.9  MAP  and  centroid  estimator  MSE  performance  as  a  function  of  average  photon  count 
per  subaperture,  K,  for  o’c/o'p  =  1/2.  Von  Karman  statistics  were  used  for  a  5  x  5 
array  of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and 
overlap  allowed. 
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Figure  5.10 


Figure  5.11 


Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
.S',  for  a  2  X  2  detector  array.  The  curves  correspond  to  a  specific  ratio  of 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 


Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
S',  for  a  4  X  4  detector  array.  The  curves  correspond  to  a  specific  ratio  of  o-^/fTp. 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 
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Figure  5.12 


Figure  5.13 


Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
jK",  for  a  8  X  8  detector  array.  The  curves  correspond  to  a  specific  ratio  of 
Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured  subapertures  with 
Lo  —  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 


MAP  and  centroid  estimator  MSE  performance  as  a  function  of  average  photon  count 
per  subaperture,  .S’,  for  a  4  x  4  detector  array.  The  curves  correspond  to  a  specific 
ratio  of  o-^/crp.  Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured 
subapertures  with  Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap  allowed. 
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5.5  Read  noise 


The  effect  of  read  noise  on  the  MSE  performance  was  examined  next.  Read  noise  was  modeled 
as  an  additive  zero  mean,  Gaussian  noise.  Read  noise  corrupts  the  detector  array  pixel  values 
resulting  in  an  increase  in  the  centroid  shift  MSE.  The  noisy  centroid  estimates  should  result  in 
MAP  shift  estimates  with  higher  MSE.  The  expected  effect  of  adding  read  noise  to  the  detection 
process  should  be  a  decreased  in  the  relative  MAP  estimator  MSE  performance. 

5.5.1  Relative  MSE  as  detector  read  noise  increases.  The  next  series  of  simulations 
examined  the  relative  MAP  estimator  MSE  performance  using  parameters  listed  in  Table  4.1  with 
read  noise,  global  tilt  removed,  overlap  allowed,  and  small  subaperture  detection  arrays.  The 
relative  MAP  estimator  MSE  performance  was  examined  for  detector  sizes  of  2  x  2  and  4x4  with 
read  noise  variances,  a\,  of  2,  5,  and  10  (photons/pixel)^.  Figures  5.14  -  5.21  show  the  results  of 
read  noise  on  the  MSE  performance  with  tilt  removed  and  overlap  allowed.  Surprisingly,  the  relative 
MAP  estimator  MSE  performance  improved  for  both  detector  array  sizes  as  the  read  noise  variance 
increased.  This  implied  that  the  MAP  estimator  was  more  robust  than  the  centroid  estimator  in 
the  presence  of  read  noise. 

As  the  read  noise  variance  increased,  there  was  more  uncertainty  in  the  MAP  estimation  at 
low  average  photon  counts  per  subaperture.  This  was  evident  by  the  random  fluxuations  of  the 
(j\  curves  in  the  relative  MSE  plots.  As  the  average  photon  count  per  subaperture  increased,  the 
effects  of  the  read  noise  were  negated  and  the  curves  smoothed  out  and  converged  towards  the 
=  0  curve  as  expected. 

5.5.2  Absolute  MSE  as  detector  read  noise  increases.  Figures  5.16,  5.19,  and  5.22  show 
the  absolute  MSE  performance  for  the  MAP  and  centroid  estimator  with  tilt  removed  wavefronts 
and  overlap  allowed.  These  results  were  not  surprising.  For  all  cases,  the  MAP  and  centroid 
estimator  MSE  performance  increased  as  the  read  noise  variance  increased.  The  MSE  performance 
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for  all  values  of  read  noise  converged  towards  the  =  0  curve  as  K  gets  large.  This  behavior 
was  expected  since  the  effects  of  read  noise  become  insignificant  as  K  gets  large.  As  the  level  of 
atmospheric  turbulence  decreased,  the  absolute  MSE  curves  also  decreased  as  expected. 


Figure  5.14  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
.S',  for  a  2  X  2  detector  array  with  read  noise.  The  individual  curves  correspond  to  a 
specific  read  noise  variance,  Von  Karman  statistics  were  used  for  a  5  x  5  array 
of  unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  al/crl  =  1,  tilt  removed, 
and  overlap  allowed. 
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Figure  5.15  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 

R",  for  a  4  X  4  detector  array  with  read  noise.  The  individual  curves  correspond  to  a 

specific  read  noise  variance,  cr^.  Von  Karman  statistics  were  used  for  a  5  x  5  array  • 

of  unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  o-^/o-p  =  1,  tilt  removed, 

and  overlap  allowed. 


Average  photon  count  per  subaperture,  K 


Figure  5.16  MAP  and  centroic^estimator  MSE  performance  as  a  function  of  average  photon  count 
per  subaperture,  .R,  for  a  4  x  4  detector  array  with  read  noise.  The  individual  curves 
correspond  to  a  specific  read  noise  variance,  (t%.  Von  Karman  statistics  were  used 
for  a  5  X  5  array  of  unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  cr^/cr^  =  1, 
tilt  removed,  and  overlap  allowed. 
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Figure  5.17  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
for  a  2  X  2  detector  array  with  read  noise.  The  individual  curves  correspond  to  a 
specific  read  noise  variance,  Von  Karman  statistics  were  used  for  a  5  x  5  array 
of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm,  aUcrp  =  1/2,  tilt  removed, 
and  overlap  allowed. 


Figure  5.18  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
if,  for  a  4  X  4  detector  array  with  read  noise.  The  individual  curves  correspond  to  a 
specific  read  noise  variance,  Von  Karman  statistics  were  used  for  a  5  x  5  array 
of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm,  (Tg/cTp  =  1/2,  tilt  removed, 
and  overlap  allowed. 
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Figure  5.19 


Figure  5.20 


MAP  and  centroid  estimator  MSE  performance  as  a  function  of  average  photon  count 
per  subaperture,  if,  for  a  4  x  4  detector  array  with  read  noise.  The  individual  curves 
correspond  to  a  specific  read  noise  variance,  Von  Karman  statistics  were  used  for 
a  5  X  5  array  of  unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  o'c/^’p  =  1/2) 
tilt  removed,  and  overlap  allowed. 


Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
if,  for  a  2  X  2  detector  array  with  read  noise.  The  individual  curves  correspond 
to  a  specific  read  noise  variance,  Von  Karman  statistics  were  used  for  a  5  x  5 
array  of  unobscured  subapertures  with  Lg  =  50m,  =  1/4,  tilt  removed,  and 

overlap  allowed. 
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Figure  5.21 


Figure  5.22 


Relative  MSB  performance  as  a  function  of  average  photon  count  per  subaperture, 
.S’,  for  a  4  X  4  detector  array  with  read  noise.  The  individual  curves  correspond 
to  a  specific  read  noise  variance,  a%.  Von  Karman  statistics  were  used  for  a  5  x  5 
array  of  unobscured  subapertures  with  Lg  =  50m,  a^/cr^  =  1/4,  tilt  removed,  and 
overlap  allowed. 


MAP  and  centroid  estimator  MSB  performance  as  a  function  of  average  photon  count 
per  subaperture,  A",  for  a  4  x  4  detector  array  with  read  noise.  The  individual  curves 
correspond  to  a  specific  read  noise  variance,  «t^.  Von  Karman  statistics  were  used  for 
a  5  X  5  array  of  unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  crUa^  =  1/4, 
tilt  removed,  and  overlap  allowed. 
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5.6  Using  incorrect  atmospheric  parameters  in  the  MAP  estimator 

Up  to  this  point,  it  was  assumed  that  the  atmospheric  parameters  were  known  with  a  high 
degree  of  certainty.  This  section  will  investigate  what  happens  to  the  MSE  performance  when 
incorrect  Fried  parameters  and  outer  scales  were  used  in  the  development  of  the  correlation  matrix 
for  the  MAP  estimator.  Since  the  Pried  parameter  has  greater  impact  on  the  mean  square  spot 
motion,  a^,  than  the  outer  scale,  more  significant  effects  should  be  observed  when  the  incorrect 

t 

Pried  parameter  is  used. 

The  graphs  in  Figs.  5.23  and  5.26  show  the  results  of  using  incorrect  Pried  parameter  (r,,) 
estimates  for  various  read  noise  variance  for  a  4  x  4  detector  array  with  tilt  removed  from  incident 
wavefronts,  and  overlap  allowed.  The  dashed  curves  represent  the  relative  MAP  estimator  MSE 
performance  when  the  correct  Pried  parameter  was  used  in  the  MAP  estimator  and  the  solid  curves 
correspond  to  the  performance  when  the  incorrect  EVied  parameter  was  used.  It  was  obvious  that 
the  choice  of  rg  has  a  significant  impact  on  the  performance  of  the  MAP  estimator.  When  rg 
was  estimated  as  higher  than  actual  (turbulence  levels  are  lower  than  estimated),  there  was  an 
improvement  in  relative  MAP  estimator  MSE  performance.  When  rg  was  estimated  lower  than 
actual,  there  was  a  decrease  in  relative  MAP  estimator  MSE  performance.  The  change  in  relative 
MSE  performance  directly  resulted  from  the  values  in  the  slope  correlation  matrix  varying  as  a 
function  of  (l/fo)®/^.  As  rg  gets  smaller,  the  correlation  values  become  larger.  As  K  becomes 
small,  the  larger  correlation  values  force  the  MAP  shift  to  zero  faster  than  if  the  correct  value  for 
rg  is  used.  Since  the  slope  correlation  matrix  determines  how  fast  the  MAP  shift  estimates  are 
driven  to  zero  as  K  gets  small,  the  choice  of  rg  will  directly  effect  the  relative  MSE  performance. 

The  graphs  in  Figs.  5.27  and  5.28  show  the  results  of  using  incorrect  outer  scale  estimates. 
Since  the  correlation  values  do  not  significantly  change  as  the  outer  scale  varied  from  10  to  100 
meters,  the  outer  scale  had  little  effect  on  the  relative  MAP  shift  estimator  MSE  performance. 
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Figure  5.23 


Figure  5.24 


Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture,  K, 
using  incorrect  Fried  parameters.  The  individual  curves  correspond  to  a  particular 
Fried  parameter,  Tq.  The  actual  Fried  parameter  was  20.58cm.  The  CCD  detector 
array  size  was  4x4.  Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured 
subapertures  with  Lq  =  50m,  d  =  9.2cm,  =  0,  tilt  removed,  and  overlap  allowed. 


MAP  and  centroid  MSE  performance  as  a  function  of  average  photon  count  per 
subaperture,  K,  using  incorrect  Pried  parameters.  The  individual  curves  correspond 
to  a  particular  Pried  parameter,  To-  The  actual  Pried  parameter  was  20.58cm.  The 
CCD  detector  array  size  was  4x4.  Von  Karman  statistics  were  used  for  a  5  x  5 
array  of  unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  cr^  =  0,  tilt  removed, 
and  overlap  allowed. 
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Figure  5.25 


Figure  5.26 


Average  photon  count  per  subaperture,  K 


Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture,  K, 
using  incorrect  Fried  parameters.  The  individual  curves  correspond  to  a  particular 
Fried  parameter,  Tq.  The  actual  Pried  parameter  was  20.58cm.  The  CCD  detector 
array  size  was  4x4.  Von  Karman  statistics  were  used  for  a  5  x  5  array  of  unobscured 
subapertures  with  Lo  =  50m,  d  =  9.2cm,  =  5,  tilt  removed,  and  overlap  allowed. 


MAP  and  centroid  MSE  performance  as  a  function  of  average  photon  count  per 
subaperture,  K,  using  incorrect  Pried  parameters.  The  individual  curves  correspond 
to  a  particular  Pried  parameter,  To-  The  actual  Pried  parameter  was  20.58cm.  The 
CCD  detector  array  size  was  4x4.  Von  Karman  statistics  were  used  for  a  5  x  5 
array  of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm,  crj^  =  5,  tilt  removed, 
and  overlap  allowed. 
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Figure  5.27  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture,  K, 
using  incorrect  outer  scales.  The  individual  curves  correspond  to  a  particular  outer 
scale  value,  Lq.  The  actual  outer  scale  was  50m.  The  CCD  detector  array  size  was 
2x2.  Read  noise  variance  was  5  (photo  events/pixel)^.  Von  Karman  statistics  were 
used  for  a  5  X  5  array  of  unobscured  subapertures  with  To  =  13.58cm,  d  =  9.2cm, 
tilt  removed,  and  overlap  allowed. 


Figure  5.28  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture,  K, 

using  incorrect  outer  scales.  The  individual  curves  correspond  to  a  particular  outer 
scale  value,  Lg.  The  actual  outer  scale  was  50m.  The  CCD  detector  array  size  was 
4x4.  Read  noise  variance  was  5  (photo  events/pixel)^.  Von  Karman  statistics  were 
used  for  a  5  X  5  array  of  unobscured  subapertures  with  Tq  =  13.58cm,  d  =  9.2cm, 
tilt  removed,  and  overlap  allowed. 
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5.7  SOR  Generation  III  detector  configuration 


The  next  set  of  simulations  examined  the  MSE  performance  when  the  subaperture  detector 
array  was  configured  to  match  those  in  the  SOR  generation  III  [19]  tests.  Each  subaperture  used 
a  4  X  4  detector  array.  In  the  modified  configuration,  the  outer  row  of  pixels  was  disabled  to 
eliminate  electronic  crosstalk  between  subapertures.  This  left  a  2  x  2  array  of  pixels  centered  on 
the  optical  axis  as  shown  in  Fig.  4.8.  The  performance  of  this  modified  array  was  compared  to  the 
performance  of  a  conventional  4x4  detector  array  at  various  levels  of  read  noise. 

Figures  5.29  -  5.31  show  that  the  standard  4x4  detector  array  had  a  better  relative  MAP 
estimator  MSE  performance  than  the  modified  detector  array  for  all  levels  of  read  noise  tested. 
The  dashed  lines  correspond  to  the  relative  MSE  performance  when  the  standard  4x4  detector 
array  was  used  and  the  solid  lines  correspond  to  the  relative  MSE  performance  when  the  modified 
detector  array  was  used.  As  the  read  noise  variance  increased,  the  relative  MSE  of  the  standard 
detector  array  outperformed  the  relative  MSE  of  the  modified  detector  array. 

However,  the  relative  MAP  estimator  MSE  performance  does  not  mean  that  the  standard 
4x4  detector  array  had  better  absolute  MSE  performance  than  the  modified  detector  array.  The 
curves  in  Fig.  5.32  show  that  the  MAP  and  centroid  estimator  MSE  performance  for  the  modified 
4x4  detector  array  was  better  than  the  standard  4x4  detector  array  for  non  zero  read  noise 
variance.  These  results  were  not  surprising  since  the  modified  detector  array  only  had  4  pixels 
effected  by  read  noise,  and  these  pixels  were  all  centrally  located  on  the  detector  array.  The 
standard  detector  array  has  16  pixels,  most  of  which  are  on  the  boundary  of  the  detector  array. 
Given  these  conditions,  the  standard  detector  array  should  have  a  higher  MSE  performance  than 
the  modified  detector  array. 
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Average  photon  count  per  subaperture,  K 


Figure  5.29  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  using  a  SOR  Generation  III  detector  with  no  read  noise.  The  pairs  of  curves 
correspond  to  a  specific  ratio  of  cr^/cTp.  Von  Karman  statistics  were  used  for  a  5  x  5 
array  of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and 
overlap  allowed. 


Average  photon  count  per  subaperture,  K 

Figure  5.30  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture, 
K,  using  a  SOR  Generation  III  detector  with  =  5.  The  pairs  of  curves  correspond 
to  a  specific  ratio  of  aHa^.  Von  Karman  statistics  were  used  for  a  5  x  5  array  of 
unobscured  subapertures  with  Lg  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap 
allowed. 
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Figure  5.31  Relative  MSE  performance  as  a  function  of  average  photon  count  per  subaperture,  K, 

using  a  SOR  Generation  III  detector  with  =  10.  The  pairs  of  curves  correspond 
to  a  specific  ratio  of  Von  Karman  statistics  were  used  for  a  5  x  5  array 

of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm,  tilt  removed,  and  overlap 
allowed. 


Figure  5.32  MSE  performance  for  the  MAP  and  the  centroid  estimator  as  a  function  of  average 
photon  count  per  subaperture,  K,  for  standard  and  modified  detector  arrays.  The 
read  noise  variance  =  10  (photons/pixel)^  and  cr^/tr^  =  1.  Von  Karman  statistics 
were  used  for  a  5  x  5  array  of  unobscured  subapertures  with  Lo  =  50m,  d  =  9.2cm, 
tilt  removed,  and  overlap  allowed. 
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VI.  Conclusions 


6.1  Introduction 

Atmospheric  turbulence  reduces  the  resolution  of  imaging  systems.  One  method  to  com¬ 
pensate  for  the  effects  of  atmospheric  turbulence  is  to  remove  the  phase  perturbations  using  an 
adaptive  optics  system.  The  two  key  components  to  remove  the  perturbations  are  the  wavefront 
sensor  and  the  deformable  mirror.  A  Hartmann  wavefront  sensor  segments  the  incident  wavefront 
over  an  array  of  subapertures  and  measures  the  centroid  shift  in  each.  Prom  this  shift  data,  the 
incident  wavefront  can  be  reconstructed  and  the  phase  perturbations  removed.  The  accuracy  of 
the  shift  measurements  will  directly  impact  the  accuracy  of  the  wavefront  reconstruction  and  the 
ability  of  the  adaptive  optics  system  to  remove  perturbations.  A  maximum  a-posteriori  (MAP) 
slope  estimator  developed  by  Sallberg  [23]  incorporated  atmospheric  statistics  and  intensity  levels 
within  the  wavefront  sensor  to  improve  centroid  shift  estimates. 

6.2  Summary  of  Methodology 

The  Hartmann  wavefront  sensor  was  modeled  in  simulation  in  order  to  evaluate  the  perfor¬ 
mance  of  the  MAP  estimator.  The  goal  was  to  evaluate  the  MAP  estimator  performance  using 
realizable  wavefront  sensor  parameters.  The  MAP  estimator  had  to  be  evaluated  in  simulation  as 
a  closed  form  solution  for  a  MAP  estimator  using  realizable  wavefront  sensor  parameters  was  not 
possible.  In  addition  to  the  wavefront  sensor  simulation,  an  implementation  for  the  MAP  estimator 
was  developed  and  phase  screen  data  was  generated  using  von  Karman  atmospheric  statistics  and 
a  Fourier  Series  phase  screen  generation  algorithm.  The  MAP  estimator  MSE  performance  was 
examined  relative  to  the  centroid  estimator  MSE  performance  as  both  the  wavefront  sensor  and 
atmospheric  parameters  were  varied. 
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6.3  MAP  estimator  performance 


The  MAP  estimator  was  found  to  perform  better  than  the  centroid  estimator  in  all  cases. 
Global  tilt  removal  had  minimal  impact  due  to  the  decreased  wavefront  correlation  statistics  result¬ 
ing  from  the  von  Karman  atmospheric  model.  Intensity  overlap  onto  adjacent  subapertures  had 
little  effect  on  the  relative  performance  of  the  MAP  estimator  as  a  result  of  the  optimized  param¬ 
eter  selection  used  in  the  wavefront  sensor  simulation.  The  relative  MSB  performance  of  the  MAP 
estimator  decreased  as  the  number  of  pixels  in  the  subaperture  detector  arrays  were  decreased. 
The  MAP  estimator  was  found  to  be  less  effected  by  read  noise  than  the  centroid  estimator.  When 
incorrect  atmospheric  parameters  were  used  in  the  development  of  the  MAP  estimator  correlation 
matrix,  it  was  found  that  incorrect  outer  scale  had  little  impact  on  the  relative  MSB  performance 
while  the  selection  of  the  Pried  parameter  significantly  affected  the  relative  MSB  performance. 
Choosing  ro  smaller  than  actual  resulted  in  increased  relative  MAP  estimate  MSB  while  choosing 
ro  larger  than  actual  resulted  in  a  decrease  in  relative  MAP  estimate  MSB.  While  the  modified 
4x4  detector  array  exhibited  better  overall  MSB  performance,  the  relative  MAP  estimator  MSB 
performance  decreased  as  expected.  In  all  cases,  the  MAP  estimator  MSB  performance  was  upper 
bounded  by  the  centroid  estimator  MSB. 

6.4  Recommendations 

The  feasibility  of  using  the  MAP  estimator  as  part  of  a  real  time  adaptive  optics  system 
should  be  studied.  Real  time  implementation  of  the  estimator  raises  several  issues.  These  include 
operation  with  a  deformable  mirror  using  a  least  squares  fitting  algorithm  in  an  adaptive  optics 
system,  the  speed  at  which  the  MAP  estimator  calculation  can  be  made,  and  the  rate  of  change  of 
the  atmospheric  parameters. 

This  thesis  examined  the  performance  of  a  MAP  estimator  using  a  ‘snapshot’  approach. 
Centroid  shift  estimates  were  made  for  individual  wavefronts  that  were  not  subject  to  a  deformable 
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mirror.  In  a  real  time  adaptive  optics  system,  the  wavefronts  will  be  altered  by  a  deformable  mirror 
prior  to  reaching  the  wavefront  sensor.  Since  the  deformable  mirror  will  significantly  change  the 
wavefront  statistics,  the  performance  of  the  MAP  estimator  should  be  examined  in  an  adaptive 
optics  system. 

Additionally,  the  MAP  estimator  inherently  requires  significant  computational  overhead.  The 
correlation  matrix  for  a  MAP  estimator  implementation  with  a  16  x  16  array  of  subapertures  will 
require  481  calculations.  As  the  atmospheric  parameters  fluctuate,  the  correlation  matrix  will  have 
to  updated.  Since  the  outer  scale  has  a  minimal  effect  on  the  correlation  values,  the  correlation  ma¬ 
trix  can  be  calculated  for  a  normalized  Pried  parameter  (t-q  =  1)  and  the  correlation  matrix  scaled 
Eis  To  changes.  Another  source  of  computational  overhead  arises  from  the  additional  512^  calcula¬ 
tions  required  to  make  the  MAP  estimate.  In  order  to  decrease  the  number  of  calculations  required 
to  make  the  MAP  shift  estimate,  approaches  using  less  correlation  data  could  be  investigated. 
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Appendix  A.  Slope  correlation  for  Kolmogorov  Statistics 


This  appendix  provides  the  derivation  for  the  equations  used  to  develop  the  slope  correlation  matrix 
using  the  Kolmogorov  structure  function.  Beginning  with  Eqn.  (3.11),  the  dot  products  x  •  dij  and 
y  ■  dij  yield  four  possible  cases  that  must  be  evaluated: 

•  When  di  and  dj  denote  x  directed  slopes. 

•  When  di  and  dj  denote  y  directed  slopes. 

•  When  di  denotes  a  y  directed  slope  and  dj  denotes  a  x  directed  slope. 

•  When  di  denotes  a  x  directed  slope  and  dj  denotes  a  y  directed  slope. 

The  first  derivation  will  be  for  the  correlation  between  subapertures  with  x  directed  slopes.  Using 
the  notation  R^x  =  Rxxixo,yo),  Eqn.  (3.11)  can  be  written  as 

^  J  J  j  j  da:da:'d2/dj^'  ^  ^S(x  +  ^)  -  S(x  -  ^)  j  rect  0)J 

^  2  ~2~  d^"  )] 

X  ^r0(0, 0)  -  ^D^(x'  -x,y'  -  y)j  .  (A.l) 

Rearranging  terms,  the  equation  can  be  written 

Rxx  =  ^  j  J  J  J  dajdx'dj/dy'rect  rect 

X  ^S{x  +  ^)S{x'  +  Six  -  ^)6ix'  +  \-Xo) 

-Six  +  \)Six'  -\-xA  +  Six  -  ^)Six'  -  ^  -  xo) 

X  ^r0(0, 0)  -  ^D^ix'  -x,y'  -  y)'j  .  (A.2) 
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Expanding  the  previous  equation: 


Rxx  =  ^  j  j  j  j  d!Eda;'dydj/'rect  Q)  rect  ^ 

X  +  ^)6{x'  +  ^-Xo)  ^r^(0,0)  -  ^D^{x'  -x,y'  -  y)^ 
-6{x  -  ^)6(x'  +  ^-Xo)  ^r^(0,0)  -  ^D4,{x'  -x,y'  -  y)^ 
-S{x  +  ^)6{x'  -^-Xo)  ^r^(0,0)  -  ^D^{x'  -x,y'  -  y)^ 
+S{x  -  ^)6{x'  -^-Xo)  (T^i0,0)  -  ^D^{x'  -x,y'  -  y)J 

The  derivation  will  make  use  of  the  Kolmogorov  phase  structure  function: 


D^{x'  -  x,y'  -  y)  =  6.88 


^{x'  -  x)^  -h  jy'  -  yif 


Using  the  sifting  property  of  dirac  delta  functions  [12]  and  recognizing  that  r,^(0, 0)  factors  out, 
Eqn.  (A.3)  can  be  written  as 


~^{hy  II  (d)  (^) 

x[((-|+x,  +  |) 


Combining  like  terms 


Rxx  ““ 


^  (^) "  (d) (^)  [2  K  +  (y'  - »)’)' 

-  ((xo  +  df  +  {y'  -  yf) «  -  ((xo  -  df  +  (?/'-  y)^)  ®]  . 
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Perform  the  following  change  of  variables: 


y'+y  j  A  I 

Sj,  =  and  Ay  =  y'-  y, 

2/'  =  Sj,  +  ^  and  y  =  Y,y-^.  (A.7) 

The  Jacobian  is  1  and  the  equation  becomes 

X  [2  {xl  +  {Ayf)  “  -  {{Xo  +  df  +  {y'  -  yf) « 

-[{xo-df  +  {y'-vff^].  (A.8) 


The  equation  can  be  simplified  to  a  single  integral  by  eliminating  one  variable.  Consider  the  terms 
in  the  previous  equation  that  depend  only  on  E„: 


tect  reel 


■‘y  ~h  2”  yo\ 


(A.9) 


Define: 


S  —  ^  1 

a  =  —■  •  and  da  =  ^dS^. 
d  d  " 


(A.10) 


The  equation  can  then  be  written  as 


d  J  da  rect  (a)  rect  ^a  — 


(A.11) 
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Since  (— Ay +2/o)/d  is  constant,  this  integral  is  the  correlation  of  the  two  rect  functions.  Referencing 
Goodman  and  using  Fourier  Transform  analysis,  the  integral  becomes 


(A.12) 


The  triangle  function  defines  the  limits  of  the  remaining  integral: 


< 


I  Ay  +  j/o|  ^ 


j/o  -  d  <  Ay  <  2/o  +  d 


otherwise 


(A.13) 


Substituting  this  into  Eqn.  (A.8) 


Rxx  — 


-  ((xo  +  d)^  +  iy'  -  yf)  “  -  {{xo  -  df  +  [y'  -  yf)  ®  ] . 


(A.14) 


Perform  the  following  change  of  variables  to  put  the  integral  in  dimensionless  quantities: 


Vo  = 
dAy  = 


d’ 

dAy 

d  ■ 


(A.15) 


The  final  equation  for  the  correlation  between  subapertures  with  x  directed  slopes  is 


Rxx  “ 


-3.44  /  d  \  3 

^  -  /  dA5(l-|Ay-y<,|) 

«  \^o/  Jyo-1 

X  [2  {xl  +  A?)  «  -  ((Xo  +  1)2  +  A|)  »  -  {{Xo  -  1)2  +  A|) 
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(A.16) 
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where  Rxx  is  a  function  of  Xo  and  yo,  the  normalized  x  and  y  separation  of  the  subapertures. 
The  development  of  the  correlation  when  both  subapertures  have  y  directed  slopes,  di,dj  =  y,  is 
identical.  The  result  can  simply  be  written  as 


_Q  44  /  \  3  /•»o+l 

X  [2  (jf2  +  A|) "  -  {{yo  +  If  +  A|) «  -  {{yo  -  1)^  +  A?) «]  .  (A.17) 

The  development  for  correlations  when  subapertures  have  slopes  in  orthogonal  directions 
follows.  The  specific  case  for  correlation  between  y  directed  slope  and  x  directed  slope  follows. 
With  di—y  and  dj=x,  the  cross  correlation  equations  can  be  found.  Using  the  notation  Ryx  = 
Rvx{xo,yo),  Eqn.  (3.11)  can  be  written  as 

^*'“  ^  / / / / ~  (d)) 

^  (i  +i~yo)-  ^(y'  -  f  -  y^))  ) 

X  ^r^(0, 0)  -  ^D^{x  -x',y-  y')j  .  (A.18) 

Rearranging  terms,  the  equation  can  be  written 

Ryx  =  ^  /  /  /  /  ^®^®^^2/d3/'rect  rect 

X  ['^(®  +  + 1  - 1/0)  -  s{x  +  ^)S{y'  -\-yo) 

-S{x  -  ^)S{y'  +  ^-yo)  +  S{x  -  ^)S{y'  -^-Vo) 

X  ^r,^(0,0)-^T>^(a;-a;',3/-y)^  .  (A.19) 
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After  substituting  in  the  structure  function  from  Eqn.(  A.4)  and  using  the  sifting  property  of  dirac 
delta  functions  [12],  the  equation  can  be  reduced  to  a  double  integration.  Unlike  the  correlation 
between  subapertures  with  slopes  in  the  same  direction,  this  equation  will  not  reduce  to  a  single 
integration. 


Ryx  " 


-3.44 


((a:'  +  +  (-|  +yo-  2/)^^  -  +  i^+yo-  yf'j 

Jj, 

-  ((»'  -  5)’  +  (-5+»«-  vf) '  +  ((»'  -  j)’  +  (f  + 11.  -  vf) 
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The  rect  functions  define  the  limits  of  integration: 


E 


yx 


-3.44 


/  1  \  3  /*a  r^o+l 

\To)  J-jJxo-i 


dx'dy 


X 


+  (-^  +  yo-  yf^  -  +  (^  +  2/o  -  y?^ 

-  ((a;'  -  +  yo-  y?^  +  ^(a;'  -  +  (^  +  j/o  -  y?'^ 


(A.21) 


Perform  the  following  change  of  variables  to  put  the  integral  in  dimensionless  quantities: 


o 

II 

and 

II 

X 

^~d 

and 

II 

-  y 

and 

dy 

(A.22) 
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The  correlation  between  y  directed  slope  and  x  directed  slope  is 


Ryx  — 


dx'dy 


(®'  +  +  (-^  +  yo-  +  i\+yo-  yf^ 

-  ((*'  -  \f  +  {-\+%-  yf^  +  -  \f  +  {\+yo-  y^^ 
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(A.23) 


where  Eyx  is  a  function  of  Xo  and  yo,  the  normalized  x  and  y  separation  of  the  subapertures.  The 
correlation  between  subapertures  with  a  x  directed  slope  and  a  y  directed  slope  is  developed  in  the 
same  manner. 


Rxy  — 


41 


(A.24) 
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Appendix  B.  Slope  correlation  for  von  Karman  Statistics 


This  appendix  provides  the  derivation  for  the  equations  used  to  develop  the  slope  correlation 
matrix  for  the  slope  correlation  matrix  using  the  Von  Karman  structure  function.  Beginning  with 
Eqn.(  3.11),  the  dot  products  x  •  dij  and  y  ■  dij  yield  four  possible  cases  that  must  be  evaluated: 

•  When  di  and  dj  denote  x  directed  slopes. 

•  When  di  and  dj  denote  y  directed  slopes. 

•  When  di  denotes  a  y  directed  slope  and  dj  denotes  a  x  directed  slope. 

•  When  di  denotes  a  x  directed  slope  and  dj  denotes  a  y  directed  slope. 

The  first  derivation  will  be  for  the  correlation  between  subapertures  with  x  directed  slopes.  Using 
the  notation  Rxx  =  Rxxixo,yo))  Eqn.  (3.11)  can  be  written  as 

Rxx  =  j  j  j  j  dxdx'dydy'  ^  (6{x  +  ^)-  6{x  -  ~)^  rect  j 

+  5  "  -  f  ‘  *»>)  (^4^)] 

X  ^r^(0,0)  -  ^D^{x'  -x,y'  -y)^  .  (B.l) 

Rearranging  terms,  the  equation  can  be  written 

Rxx  =  ^  J  f  j  J  ^^dx'dydy'iect  rect  ^  ^ 

X  +  ^)6{x'  +  |  -  x<,)  -  «(a;  -  ^)6{x'  +  ^-Xo) 

-6{x  +  |)<5(x'  -^-Xo)  +  Six  -  ^)S{x'  -  I  -  aJo) 

X  ^r0(0, 0)  -  ^D^ix'  -x,y'  -  y)j  .  (B.2) 
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Expanding  the  previous  equation, 


R 


XX 


X  [<5(3^  +  +  \-Xo)  ^r^(0,0)  -  ^D^{x'  -x,y'-  y)j 

-S{x  -  ^)6{x'  +  ^-Xo)  ^r^(0,0)  -  ^D4,{x'  -x,y'  -  y)^ 
-Six  +  ^)Six'  -^-Xo)  (r^iO,  0)  -  ^D^ix'  -x,y'  -  y)^ 
+Six  -  ^)Six'  -^-Xo)  ^r^(0,0)  -  ^D^ix'  -x,y'  -  y)^ 


(B.3) 


The  derivation  will  make  use  of  the  von  Karman  phase  structure  function  [30]: 


D^ix)  = 


3.089  6 

(27r)3  5 


(B.4) 


where  r[-]  is  the  Gamma  function,  li’s/el-]  is  a  modified  Bessel  function  of  the  second  kind  of  order 
5/6,  Lo  is  the  outer  scale  parameter,  Tq  is  the  Pried  parameter,  and  |f|  =  y/{x'  -  x)^  +  {y'  -  y)^. 

Using  the  sifting  property  of  dirac  delta  functions  [12]  and  recognizing  that  r^(0,0)  factors 
out,  Eqn.  (B.3)  can  be  written  as 


Rx 


=  0.08663  (^)  •  1  /  /  drt'iect  (^)  rect  (S^) 


,  r[l/61  /'V(i,-dP+(g'-tf' 


1 

TTS 


K, 


5/6 


27r 


y/jxo-d)^  +  {y'  -yy 


+ 


r[i/6i  /  ^/(^7W+F^'\  *  „ 

L.  )  ’‘‘I’ 


27r 


^/{xo  +  d)‘^  +  jy'  -yY 


_2  (l  nm  fy/xl  +  jy'-yf 
,  TTO  V  Lo 


K, 


6/6 


27r 


Vxl  +  jy'  -  yf 


(B.5) 
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Perform  the  following  change  of  variables: 


V  y'+y  j  A  / 

Si,  =  and  Ay  =  y  -y, 

'/  =  '^y  +  ^  and  y  =  'Ey-^. 


(B.6) 


The  Jacobian  is  1  and  the  equation  becomes 


Rxx  =  0.08663 


iS'ill 


/  E  -  ^ 
dA„dEi,rect  I  ^  ^  ^ 


rect 


'^y  +  ^-yo' 


r[i/6]  f  ^/xl  +  jAy)^' 


7ra 


Lo 


K, 


5/6 


27r 


V^l  +  jAyr 


r[l/6]  f^ixo  +  dy  +  {Ayr' 


i 

Tra 


K. 


5/6 


27r 


r[l/6]  f  ,/iXo-d)^  +  (Ay  f 


1 

Tra 


K> 


5/6 


27r 


^iXo+d)^  +  {AyY' 

Lo 

^{Xo  -  dy  +  {AyY 


(B.7) 


The  previous  equation  can  be  simplified  to  a  single  integral  by  eliminating  one  variable.  Consider 
the  terms  that  depend  only  on  E„: 


/ 


dS^irect 


'v 

2 


rect 


-  i/o 


(B.8) 


Define: 


a  = 


V  _  Am. 

_  ^y  2 


and  da  =  -dE„. 
d 


(B.9) 


Equation  (B.8)  can  then  be  written  as 


d  J  darect  (a)  rect 


(B.IO) 
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Since  (— Ay  4-yo)/<^  is  constant,  this  integral  is  the  correlation  of  the  two  rect  functions.  Referencing 
Goodman  and  using  Fourier  Transform  analysis,  the  integral  becomes: 


/ds:. 


(B.ll) 


The  triangle  function  defines  the  limits  of  the  remaining  integral: 


“(^)  = 


(■^.-b.Av+yo[)  y^-d<^^y<yo  +  d 

0  otherwise. 


(B.12) 


Substituting  this  into  Eqn.  (B.7): 


Rx 


/  r  \  3  1  rvo+d  / 


-  1  -  Ay  +1 


— L. —  2* — r, — 


Tra 


1 

Tra 


fViXo  +  dr  +  iAyrW^ 

i.  j 

Lo 

/V(a:.-d)2  +  (Ay)2 

V^(a:o  -  d)2  +  (Ay)2 
Lo 

(B.13) 


Perform  the  following  change  of  variables  to  put  the  limits  of  integration  in  dimensionless  quantities: 


and 


Vo 


and 


dAy  = 


dAy 

d  ' 


(B.14) 
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The  final  equation  for  the  correlation  between  subapertures  with  x  directed  slopes  is: 


■  -©'ira/::**-"-''-"’ 


\  Lo  )  Lo 


dV^7+W+T^Y  ^  ^/Vixo  +  ir  +  jAy^ 
Loj  Lo 

dy/{Xo-l?  +  iAy)^\  '  L_M^o-l)2  +  (Aj;)2 

- Z -  - 


(B.15) 


where  Rxx  is  a  function  of  Xq  and  yo,  the  normalized  hatx  and  y  separation  of  the  subapertures. 
The  development  of  the  correlation  between  subapertures  with  y  directed  slopes,  di  and  dj  =  y,  is 
identical.  The  result  can  simply  be  written  as 


=  0.08663  (^)  *  i  (£^)  £"■  dA,(l  -  |A.  -  »,|) 


\  Lo  )  Lo 


V(yo  +  l)^  +  (Ax)n "  LdV{yo  +  ir  +  {A,r 


dy/iyo-ir  +  iAxr] ^_d^{yo-ir  +  {Axr 

- F -  ^6/6  27r - - - 


(B.16) 
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The  correlation  between  subapertures  with  y  directed  slope  and  x  directed  slope,  di=y  and 
dj=x,  follows.  Using  the  notation  Ryx  =  RyxixoiVo)-  Eqn.  (3.11)  can  be  written  as 

Ryx  =  J  j  J  j  dxdx'dydy' (s{x  +  ^)-S{x-^)^Tect(^'^^ 

{i  +!-»“)-  ■*(»'  -  5  -  >'■>))  ) 

X  (r,(0,0)  -  iD,(i  -x’,y-  9'))  .  (B.17) 


Rearranging  terms,  the  equation  can  be  written 

Ryx  =  ji  J  J  j  J  da;da;'d2/d2/'rect  rect 

X  [<5(2:  +  ^)Siy'  +  ^-yo)-  i(x  +  |)%'  -\-yo) 

-6{x  - + 1  -  j/o)  +  Kx  - 1)%'  -\-yo) 

X  (B.18) 

Use  the  von  Karman  structure  function  defined  in  Eqn.  (B.4)  and  the  sifting  property  of  delta 
functions 
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The  rect  functions  define  the  limits  of  integration: 


R 


■yx 


=  T‘LU 


Tre 

r[l/6]  / 

1 

.  t 

a/(®'-4)H(4+vo-vP  \  ® 

TTB  \ 

) 

r[i/6]  / 

\/(x’+^)^+{-i+yo-yy  \ 

tts  \ 

r[i/6] , 

J 

.  1 

'  ^/ix'+i)H(i+yo-y)^  \  ” 

\2  \  6 


^5/6 


Hi+Vo-y)^ 


•^6/6 

^5/6 


2^\/(.^'+i)^+(i+yo-^ 


(B.20) 


Perform  the  following  change  of  variables  to  put  the  limits  of  integration  in  dimensionless  quantities: 


Xo  A  ~  Vo 
Xo  =  -T  and  Vo  =  -j, 


x'  =  ^  and  di'  =  ^ 
a  a 

y'  =  2 


(B.21) 


Factoring  out  common  terms,  the  correlation  between  subapertures  with  y  directed  slope  and  x 
directed  slope  is 


Ryx  — 


+i/o-g)° 

+yi,-y)^ 


(B.22) 
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where  Ryx  is  a  function  of  Xq  and  j/o,  the  normalized  hatx  and  y  separation  of  the  subapertures. 
The  correlation  between  subapertures  with  x  directed  slope  and  y  directed  slope  is  developed  in 
the  same  manner. 
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Appendix  C.  Shot  limited  performance 

This  appendix  derives  a  centroid  shift  estimator  optimized  for  a  quadcell  detector  with  shot  noise. 
Equations  developed  by  Tyler  and  Fried  [27]  were  used  for  this  derivation.  Given  incident  plane 

Y 

a 

1  2 

^ - ►X 

3  4 

Figure  C.l  Quad  cell  detector  configuration 

waves  in  the  subaperture  and  no  read  noise  in  the  detection  process,  the  authors  derived  a  shift 
estimator  that  was  optimized  for  shot  noise.  A  single  lens  imaging  system  with  a  detector  located 
at  the  focal  length  forms  an  intensity  distribution  centered  on  the  optical  axis  when  when  incident 
wavefront  is  planer.  Shot  noise  in  the  detection  process  will  result  in  centroids  shifted  off  the 
optical  axis.  The  distance  the  centroid  deviated  from  the  optical  axis  is  the  shot  induced  error.  If 
the  centroid  shift  errors  are  assumed  to  be  small,  a  linear  relationship  can  be  found  between  the 
detected  centroid  position,  x,  and  the  difference  in  intensities  detected  on  each  half  of  the  detector 
in  the  x  and  y  directions.  A/®  and  Aly.  Since  the  actual  centroid  distributions  are  centered  on 
the  optical  axis,  the  equation  developed  in  this  appendix  will  be  the  shot  limited  estimate  for  the 
centroid  shift. 

The  following  derivation  is  for  the  shift  in  the  x  direction.  A4  is  the  sum  of  the  detected 
signal  on  quadrants  2  and  4  minus  the  sum  of  the  detected  signal  on  quadrants  1  and  3. 

The  expected  intensity  difference  is  defined: 

|»00 

A4=/  d4F(4  +  4j- 

Jo 


i: 


daF(4  +  4J, 


(C.l) 
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where  F  was  the  one  dimensional  point  spread  function  (PSF): 

Fifx)  =  j  dxi{x,0)Tix,0)exi>^^^f‘\  (C.2) 

First  the  Fourier  transform  of  the  irradiance  on  the  detector,  I,  is  defined: 
i{x,y)  =  J  yd/,d/t,C»,(/„/„)exp-2-JW*+«'/v),  (C.3) 

where  Ooifx,fy)  is  the  object  irradiance.  Using  a  point  source  irradiance,  Ooifx,fy)  is  defined: 

Ooifxj  fy)  =  fy)}  (C.4) 

dyj 

where  lo  is  the  total  intensity  on  the  detector  and  dw  is  the  width  of  the  wavefront  sensor  subaper¬ 
ture.  Solving  for  7  in  Eqn  (C.3): 

=  (C.5) 

dxu 

Next,  the  optical  transfer  function  T{fa,  fp)  must  be  approximated.  The  result  is  simply  the 
autocorrelation  of  the  pupil  function: 


nx,y)  =  J  y'dAd/j,|P(^,/j,)pexp-2-JW.+y/v) 

=  :F-^[Pif.,fy)]*:F-^[PifxJy)] 

=  P{x,y)-kP{x,y),  (C.6) 
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where  T  is  the  Fourier  transform  and  the  pupil  function  for  a  square  aperture  is  defined: 

P{x,y)  =  rect^^)rect(J^),  (C.7) 

where  the  rect  function  was  defined  in  Eqn.  (2.5).  The  optical  transfer  function  in  Eqn.  (C.6)  can 
be  written  as: 

-  ^)  (l  -  ^)  for  la:|<da,  and  \y\<dw.  (C.8) 

Solving  for  F  in  Eqn.  (C.2),  the  1  dimensional  PSF  became: 

F(/^)  =  I dxIo{l-  ^)exp2’^^“/« 

—  sine  {dyjfx)' 

Since  the  centroid  displacement  was  small,  the  authors  used  the  first  order  term  of  a  Taylor 
series  expansion  about  Xo  for  Eqn.  (C.l)  in  order  to  relate  the  difference  in  intensity  AJj,  to  the 
object  displacement  Xg. 

A4  =  -2gp’(0).  (C.IO) 

Solving  for  F{0)  in  Eqn.  (C.9): 

F(0)  =  Igdg,.  (C.ll) 
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The  final  form  of  the  centroid  estimator  can  now  be  found.  Solving  Eqn.  (C.IO)  for  the 


centroid  shift  in  the  x  direction,  Xq' 


-A4  Xfi 
lo  ^dw  ’ 


(C.12) 


where  z  is  the  focal  length  of  the  system  fi,  A  is  the  average  wavelength,  is  the  width  of  the 
wavefront  sensor  subaperture,  AIx  is  the  the  sum  of  the  detected  signal  on  quadrants  2  and  4 
minus  the  sum  of  the  detected  signal  on  quadrants  1  and  3,  and  lo  is  the  total  detected  signal.  The 
centroid  shift  in  the  y  direction  is  found  in  the  same  manner. 


Vo  = 


Iq  2dt(i 


(C.13) 


where  Aly  is  the  the  sum  of  the  detected  signal  on  quadrants  3  and  4  minus  the  sum  of  the  detected 
signal  on  quadrants  1  and  2.  The  centroid  shift  estimate  can  then  be  written  as 


^  =  ajox  +  Voy. 


(C.14) 


C-4 


S?  S?  S?  S?  SS  S?  S?  S? 


Appendix  D.  Wavefront  sensor  source  code  listing 


This  Appendix  contains  the  source  code  listing  for  the  wavefront  sensor  and  the  functions  it  calls. 


D.l  Wavefront  sensor  code 

^HiH:^^:**S|:*^^:^:^:^HH^*^t^l******************************** 

Troy  B  Van  Caeter  GE-97D 
Hartmann  Wavefront  Sensor  Simulation 
Last  Modified:  SO  Jun  1997 

iliit:lti^:if}lf*iH::H:S‘**<lt************************S‘********* 

This  function  returns  a  vector  xest  with  x  and  y  shift  estimates  and 

a  vector  xact  with  actual  x  and  y  shift  values  for  the  centroid  10 

location  in  each  subaperture.  The  function  is  called  as  follows: 

[xest, xact]  =  wfs(screenl, naps, foe, diam, lambda, ccdnumb, 
m,  K,  tilton,  dutyc,  solap ) 

function  [xest, xact]  =  wfsl(screenl,  . . . 

naps,  . . . 
foe,  . . . 
diam,  . . . 

lambda,  ...  20 

ccdnumb,  . . . 
rn,  ...  . 

K,  ... 
tilton,  . . . 
dutyc,  ... 
solap) 

%  screenl  is  the  phase  screen  in  the  aperture  of  the 
%  wavefront  sensor. 

%  30 

%  Parameters  in  the  function  call 

%  naps  -  N  where  N*N  is  the  total  number  of  subapertures 

%  foe  -  subaperture  focal  length 

%  diam  -  subpaerture  pupil  diameter 

%  lambda  -  average  wavelength 

%  ccdnumb  -  M  where  M*M  is  the  number  of  detector  pixels  per  subaperture 

%  m  -  read  noise  variance 

%  K  -  average  photon  count  per  subaperture 

%  tilton  -  one=remove  tilt  zero=retain  tilt 

%  dutyc  -  duty  cycle  of  the  detector  40 

%  solap  -  subaperture  overlap  one— allowed  zero=not  allowed 

% 

%  Variables  used  in  this  function 

%  nft  -  the  size  of  the  matrix  containing  the  intensity  distribution 
%  data  for  one  subaperture  detector 

%  Na  -  the  array  size  of  the  phase  screen  in  each  subaperture 
%  nn  -  the  overall  array  size 

%  amattot  -  the  matrix  containing  complex  amplitude  data  for  all 
%  subapertures  in  the  wavefront  sensor  array 

%  mattot  -  the  matrix  containing  the  intensity  distribution  for  all  50 

%  subapertures  in  the  wavefront  sensor  array 

%  sapmask  -  subaperture  mask  for  intensity  overlap 
%  matl,matZ  -  used  to  find  the  estimated  centroid  position 
%  sic  -  used  to  make  the  intensity  pattern  an  even  multiple 
%  of  the  CCD  detector  array 

%  delx,dely  -  vectors  used  in  the  moment  calculation 
%  dphix,dphiy  -  used  in  the  actual  centroid  calculation 

%  *:K*i|<*>l<**H<N<i|i**********»****>l<*>|iiKi|i***i|<****N<*****1<********* 

%  function  begins  60 
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%  calculate  some  preliminary  items 

Na=15;  %This  is  the  number  of  samples  in  the  subaperture 

%  Determine  the  number  of  pixels  in  the  transform  matrix  window  that  are  70 

%  required  to  physically  occupy  the  space  on  the  detector  array  in  each 
%  subaperture. 

nft=floor((diajn*64)/((Na— l)*lambda*(foc/diam)))+l; 

%  Make  sure  nft<=63,  otherwise  the  indexes  to  matricies  are  all  messed  up. 

%  If  nft  is  too  large,  increase  the  subaperture  “window”  size.  If  nft 
%  cannot  be  made  small  enough  by  increasing  Na,  the  function  will  crash. 

while  nft>63  &  Na<40  80 

Na=Na+l; 

nft=floor((diam*64)/((Na— l)*lambda*(foc/diam)))+l; 
end 

%  Force  the  size  of  the  transform  matrix  window  to  be  odd  and  center  it  on 
%  pixel  (33,33).  This  places  the  maximum  intensity  value  at  the  center  of 
%  each  subaperture  and  ensures  the  fftshift  command  will  properly  orient  the 
%  maximum  to  the  upper  left  comer  before  using  the  fftS  command. 

if  floor(nft/2)==nft/2  90 

nft=nft— 1; 
end 

%  Create  a  mask  for  the  subaperture.  If  no  overlap  is  desired,  mask  the 
%  ovelap  area  with  zeros. 

if  solap==l 

8apmask=ones(63) ; 
else 

8apmask=zeros(63);  100 

sapmask(32— (nft-l)/2:32+(nft— l)/2,32-(nft-l)/2:32+(nft-l)/2)=ones(nft); 
end 

%  Scale  the  photon  count  per  subaperture  by  the  duty  cycle  of  the  detector  array. 

K=K*dutyc; 

%  Remove  tilt  from  phase  screen  if  required. 

If  tilton==l  110 

screen  1 =tiltofF(screen  1 ) ; 
end 

%  Initialize  the  detector  array  (amattot)  to  zero. 

amattot=zeros(nft*naps+(64— 1— nft)); 

%  Increase  the  number  of  samples  in  the  phase  screen  so  there  are  Na*Na 
%  samples  for  each  subaperture  screenl=resize(screenl,naps*Na). 

120 

%  Hi****************************************************** 

%  compute  complex  amplitude  on  the  detector 

% 

p=l; 

for  hc=l:naps 
for  vc=l:naps 

%  Window  the  phase  screen  data  into  the  64*64  array. 

130 

matl=zeros(64); 

matl(33-(Na-l)/2:33+(Na-l)/2,33-(Na-l)/2:33+(Na-l)/2)=  . . . 
exp(screenl((hc— l)*Na+l:(hc-l)*Na+Na,  ... 

(vc-l)*Na+l:(vc— l)*Na+Na).*(l*sqrt(— 1))); 

%  Calculate  the  actual  x  and  y  positions. 

dphix=(sum(screenl((hc-l)*Na+l:(hc— l)#Na+Na,(vc— l)*Na+Na))-. .. 

sum(screenl((hc— l)*Na+l:(hc-l)*Na+Na,(vc— l)*Na+l)))/Na; 
dphiy=(sum(screenl((hc-l)*Na+Na,(vc-l)*Na+l:(vc— l)*Na+Na))-. . .  140 

Bum(screenl((hc-l)*Na+l,(vc-l)*Na+l:(vc— l)*Na+Na)))/Na; 
xact(p:p+l,l)= [dphix*foc*lambda/  (2>i‘pi*diam)  ;dphiy*foc*lambda/  (2'i<pi*diam)] ; 
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P=P+2; 

%  Take  FFT  of  the  wavefront  in  the  lubaperture  to  get  complex  amplitude  data. 

inat2=fltshift(fft2(fftshift(matl))); 

%  Add  complex  amplitudes  to  amattot  matrix. 

150 

amattot((hc— l)*nft+l:(hc— l)*nft+64— l,(vc— l)*nft+l:(vc— l)*nft+64— 1)=. . . 
amattot((hc— l)*nft+l:(hc— l)*nft+64— l,(vc— l)*nffc+l:(vc— l)*nft+64— 1)+. . . 
mat2(2:64,2:64).*sapmask; 

end 

end 

%  **sne*>Hf^t>ifiHf^ctHfSHfiniHtilHfilttli:iH‘****:**iHit*******>lf*<Klr*********** 

%  Calculate  intensity  on  the  detector  array 

mattot=amattot.*coiy(amattot); 

%  Determine  detected  intensities  for  each  subaperture 
%  and  find  the  estimated  x  and  y  offsets 

%  scl  is  an  integer  multiple  of  the  ccd  array  size.  The  intensity  data  in 

%  each  subaperture  is  resized  from  nft*nft  to  the  smallest  integer  multiple  170 

%  of  ccdnumb  larger  than  nft. 

scl=2; 

while  scl*ccdnumb  <  nft 
scl=scl+l; 

end 

for  hc=l:napa 
for  vc=l:napa 

180 

%  Pull  off  individual  subaperture  intensity  patterns  from  mattot. 
matl=2eros(nft); 

matl=mattot((64— nft-l)/2+(hc-l)*nft+l:(64— nft— l)/2+(hc-l)*nft+nft,. . . 

(64— nft— l)/2+(vc— l)*nft+l:(64— nft— l)/2-f(vc— l)*nft4-nft); 
mat2=matl./niax(max(matl));  %data  in  mat2 

%  Find  the  detected  photon  count  for  each  pixel  in  the  subaperture  detector 
%  array. 

190 

matl=resi2e(mat2,scl*ccdnumb);  %data  in  matl 
mat2=zeros(ccdnumb); 

for  x=l:ccdnumb 
for  y=l:ccdnumb 

mat2(x,y)=sum(sum(matl((x— l)*scl+l:(x— l)»scH-8cl,. . . 

(y-l)*scl+l:(y-l)*scl+scl))); 

end 

end  %data  in  matS 

200 

%  Normalize  the  subaperture  photon  count  to  K. 
total=sum(sum(mat2)); 

matl=((mat2.*K)./total);  %data  in  matl 

%  Account  for  shot  noise. 
mat2=:poisson2(matl); 

matl=mat2;  %data  in  matl 

210 

%  Add  read  noise. 

mat2=matl+rn*randn(ccdnumb);  %data  in  matS 

%  Calculate  x  and  y  offsets  using  a  moment  calculation. 

delx=(((l:ccdnumb)— (ccdnumb+l)/2)*diajn/ccdnumb) ' ; 
if  sum(sum(mat2))==0 
xoff(hc,vc)=0; 

yoff(hc,vc)=0;  220 

else 
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xoff(hc,vc)=8um(mat2)*delx/sum(8um(mat2)); 
yoff(hc,vc)=8um(mat2 ' )  ♦delx/sua(suin(mat2) ) ; 
end 

end 

end 

%  4ii|<«H<*N<*4<**>K*H<*4<**4<*4<>K*4<>|i***>|i******>I>4<**i|ii|>M<4<*4'«*«**i|»I>**** 

%  Stack  shift  estimate  results  into  a  column  vector  230 

p=l; 

for  hc=l:naps 
for  vc=l:naps 

xest(p:p+l,l)=[xoff(hc,vc);yofif(hc,vc)]; 

P=P+2; 

end 

end 

%  end  of  function  240 

%  ^^sltili**:):************************************************ 


D.2  Poisson  number  generator 

%  This  program  takes  an  n*n  matrix  and  returns 
%  an  n*n  matrix  whose  values  are  a  random  deviate 
%  drawn  from  a  poisson  distribution  generated  from 
%  the  value  of  each  entry  using  rand 
%  as  a  source  for  uniformly  distributed  deviates 
%  This  program  was  generated  from  Numerical  Recipe 
%  The  Art  of  Scientific  Computing  1986,  p.207. 

function  [newmat]  =  poisson2(mat) 

10 

rand( '  seed '  ,8um(100Kclock)); 
n=leng^h(mat); 
for  x=l:n 
for  y=l:n 
pmean=mat(x,y); 

if  pmean<12  %use  direct  method 

gpois=exp(— pmean); 
empoi8=— 1; 
tpois=l; 

empois=empois+l;  20 

tpois=tpois*rand; 
while  tpois>gpois 
empois=empois+l| 
tpois=tpois*rand; 
end 

el8e  %«»e  the  rejection  method!^ 

sqpois=sqrt(2*pmean); 
alxmpois=log(pmeaii); 

gpois=pmean*alxmpoi8— gammaln(pmean+l); 

ypoi8=tan(pi*rand);  %Y  is  a  deviate  from  a  Lorentzian  comparison  30 

empois=sqpoi8*ypoi8+pmean; 
while  empois<0  %  continue  until  empois>=0 

ypoi8=tan(pi*rand);  %y  t»  o  deviate  from  a  Lorentzian  comparison 
empoi8=sqpoi8*ypois+pmean| 
end 

empoi8=floor(empois); 

tpois=0.9*(l+ypois*2)>tiexp(empois*alxmpoi8— gammaln(empois+l)— gpois); 
while  rand  >  tpois 
ypois=tan(pi4^rand); 

empois=sqpois*ypoi8+pmean;  40 

while  empois<0 
ypois=tan(pi*rand); 
empois=sqpois*ypois+pmean; 

end 

enipois=floor(empois); 

tpois=0.9*(l+ypois"2)*exp(empois*alxmpoi8— ganimaln(empois+l)— gpois); 
end 
end 

newniat(x,y)=empois; 
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50 


end 

end 


D.3  Phase  screen  tilt  removal 


%  **  Troy  B  Van  Caster  GE-97D 
%  ♦♦  Phase  screen  tilt  removal 
%  ♦*  Last  modified  15  May  97 

^  :fc:|e3(c:^)|(3)c:|e:fc]|t:fc:|e:|c:^)|c:|e:f3|c3f)|t:)()|c:tt%]t<itc:tc9tc])e:)c^]t<i4c 

%  This  subroutine  removes  the  tilt  from  a  phase  screen.  10 

function[tilto]  =:  tiltoff(tiH) 

l=Iength(tilt); 

cent=(l+l)/2| 

sumx=0; 

sumy=0; 

sumxx=0; 

sumyy=:0;  20 

for  x=:l:l 
for  y=l:l 

sumx=sumx+tilt(x,y)*(x— cent)/cent; 
sumy =sumy +tilt  (x,y )  »(y — cent )  /cent ; 
sumxx=sumxx+((x— cent)/cent)*2; 
sumyy=sumyy+((y— cent)/cent)*2; 

end 

end 

30 

tx=sumx./sumxx; 

ty=sumy./sumyy; 

for  x=l:l 
for  y=l:l 

tilto(x,y)=tilt(x,y)-tx*(x-cent)/cent-ty*(y-cent)/cent; 

end 

end 
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Appendix  E.  Correlation  generation  source  code  listing 


This  Appendix  contains  the  source  code  listing  the  correlation  matrices. 


E.l  Source  code  for  building  Kolmogorov  correlation  matrices 


%  **  Troy  B  Van  Caster  GE-97D  ♦* 

%  **  Kolmogorov  slope  correlation  ** 

%  *♦  matrix  generator  ** 

%  Last  Modified  S  Jun  97  ** 

%  4<*4<+*>l<>l<**4<!|<4»l«l<4<*>l<«*>l<>l>«*************** 

% 

%  This  builds  a  slope  covariance  matrix  using  the  kolmogorov 
%  structure  function.  The  program  calculates  the  first  row 
%  values  then  uses  symmetry  to  fill  in  the  rest  of  the  matrix. 

%  The  program  builds  x  offset  and  y  offset  matricies  as  well 
%  as  a  ref  matrix,  xoff  and  yoff  are  subap*subap  in  size  and 
%  contain  normalized  x  and  y  offsets  from  xl,yl  to  x2,y2.  The 
%  ref  matrix  is  used  to  fill  in  the  rest  of  the  correlation 
%  matrix.  It  is  indexed  using  the  x  and  y  offsets  and  contains 
%  a  flag  to  the  appropriate  first  row  element  to  place  in  the 
%  matrix.  Called  by  R  =  kcorrmat(asize,d,ro) 

%  This  program  calls  the  following  functions: 

%  kolxx  -  for  correlation  of  slopes  that  are  in  the  same  direction 
%  for  XX  slopes;  val=quad8(’kolxx’,-l+xoff,l+xoff, [],[], xoff, yoff) 

%  xoff  is  the  X  offset  normalized  by  d 

%  yoff  is  the  y  offset  normalized  by  d 

%  limits  of  integration  are  from  -1  to  1 

% 

%  kolxyfun  •  returns  xy  slope  correlation 
%  calls  function  intgrate.m  and  kolxy.m 

%  xoff  is  the  X  offset  normalized  by  d 

%  yoff  is  the  y  ^set  normalized  by  d 

%  limits  of  integration  are  from  -If  2  to  1)2 
% 

%  kolxy  -  determines  xy  slope  correlation 
%  val=quad8(’kolxy’,xoff-i/2,xoff+l/2,f ],[ ],y, xoff, yoff) 

%  must  be  used  in  a  2D  integration  format 

%  ****+♦*♦*♦*+*♦*♦♦+♦*♦♦****+*+♦***♦♦***♦♦***♦*♦*♦♦*♦*♦♦♦*♦*♦***♦*** 

function  R  =  kcorrmat(asize,d,ro) 
subap=asize*2;  %  length  of  covariance  matrix 

%  build  X  and  y  offset  matricies  and  the  reference  matrix 

index=0; 

for  x=l;asize 

for  y=l:asize 

index=index+l; 

ref  (y  ,x) =mdex*2 — 1 ; 

xtemp(x,y)=(y-x); 

end 

end 

for  x=l:aaize 
for  y=l:asize 

xoff(aaize*x— asize+l:asize*x,asize*y— asize+l:asize*y)=xtemp; 

end 

endM 

for  x=l:asize‘2 
for  y=l;asize*2 

yofr(x,y)=(floor((y-l)/asize)-floor((x-l)/asize)); 
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end 

end 


*s|<*>l<*>K**>K4<**>l<**»*>K****4<>K4<*>K4i>K*>K***4>*4<>l<4«K4i***>K*>K*4<***N<***H<*>K**H>**K>l< 

%  calculate  the  first  row  of  the  correlation  matrix 

%  Si:^f^!^^:^^S‘********S‘**S‘>i‘**S(*S:*****St**^*StStStStSiS^!S^’¥StittilfiiSSt**********>S*****>lf 

x=l; 

for  y=l:subap*2 

if  (floor(y/2))  "=  (y/2)  %its  an  xx/yy  correlation 

R(x,y)=quad8('kolxx',yoff(l+floor(x/2),l+floor(y/2))-l,yofr(l+floor(x/2),l+floor(y/2))+l,. . . 

[],[], xoff(l+floor(x/2),l+fioor(y/2)),yoff(l+floor(x/2),l+floor(y/2))); 
else  %its  an  xy/yx  correlation 

%check  to  see  if  the  subapertrues  are  in  the  same  row/ column 
if  xofF(l+floor(x/2),floor(y/2))*yoff(l+floor(x/2),floor(y/2))==0 
R(x,y)=0; 
else 

R(x,y)=kolxyfun(xoff(l+floor(x/2),floor(y/2)),yoff(l+floor(x/2),floor(y/2))); 

end 

end 

end 

%  Fill  in  the  rest  of  the  upper  diaganol  matrix 
%  start  with  y  rows 

for  x=2:2:subap*2 
for  y=x:subap*2 

if  (floor(y/2))  ==  (y/2)  %its  an  xx/yy  correlation 

R(x,y)=R(l,  ref(abs(yoff(floor(x/2),floor(y/2)))4-l,.  •• 

abs(xoff(floor(x/2),floor(y/2)))+l)); 
else  yits  an  xy/yx  correlation 

R(x,y)=R(l,ref(ab8(yoff(floor(x/2),floor(y/2)+l))+l,. . . 

abs(xoff(floor(x/2),floor(y/2)+l))+l  )+l); 

%check  the  x  and  y  offsets  to  determine  sign  of  cross  terms 
%remember,  they  are  all  negative  to  start  with 
%if  the  xoff  or  the  yoff  is  negative,  make  the  cross  term  positive 
if  (  xoff(floor(x/2),floor(y/2)+l)  )»(  yofT(floor(x/2),floor(y/2)+l)  )  <  0 
R(x,y)=-l*R(x,y); 
end 
end 
end 
end 

%  add  X  rows 
for  x=3:2:subap*2 
for  y=x:8ubap*2 

if  (floor(y/2))  (y/2)  %its  an  xx/yy  correlation 

R(x,y)=R(l,  ref(ab8(xoff(floor(x/2)+l,floor(y/2)+l))+l,. . . 

abs(yoff(floor(x/2)+l,floor(y/2)+l))+l  )); 
else  %its  an  xy/yx  correlation 

R(x,y)=R(l,ref(ab8(xoff(floor(x/2)+l,floor(y/2)))+l,. . . 

abs(yoff(floor(x/2)+l,floor(y/2)))+l)+l); 

%  check  the  x  offset  and  if  negative,  make  the  correlation  value  negative 
if  (  xoiF(floor(x/2)+l,floor(y/2))  )*(  yoff(floor(x/2)+l,floor(y/2))  )<  0 
R(x,y)=-l*R(x,y); 
end 
end 
end 
end 
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%**************^H,*****S‘*****************************S:************** 

%mirror  the  matrix 
R=R+triu(R,l) ' ; 

R=R»(l/d)*2*(d/ro)-(S/3); 

%«***K»K>K***>K4<**iKiK*>l<«i|<4<>K*i|<*>K******sf»*«**************i|c«*******i|i;«i<i**>l<*  ■ 

_  130 


%  **  Troy  B  Van  Caster  GE-97D  ♦* 
%  *♦  Kolmogorov  slope  correlation  ♦* 

%  *♦  matrix  generator  *♦ 
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%  Last  Modified  5  Jun  97  ** 

%  ****♦***♦*****♦********♦***♦*♦»***♦** 

%  He *****>(<** 

10 

function  val=kolxx(y,xoff,yoff) 

val=3.44*(l— abs(yoflF— y)).*(. . . 

((xoff+l)-2+y.-2).-(5/6)+... 

((xoff-l)-2+y.-2).*(5/6)-... 

2*(xoflF‘2+y."2).‘(5/6)  ); 


%  Hci|eHeHe)(eH(>l<H(HcH(>(ci|tH(i(eH(>(e4:>|e>(eete>(e>)e>l<i(e>l<ifeel<ittH(HeHe4e)((e)<i((*il< 

%  HeHeHe>KHeHei|eHeHe*HeHei|eiKH<>Kel'>KMe>leeK*4eHe**HeHei|e4e>|e*e|ei|e*eie* 

%  »*  Troy  B  Van  Caster  GE-97D  *♦ 

%  **  ZD  integrator  for  xy  and  yx  ** 

%  **  slope  correlation  terms  ** 

%  eieHe  Lost  Modified  3  May  97  ** 

%  HeHeHeHeHeHee|eHei|eHei|eHei|eHeHei|e4ei|eHee|ei|eHeHeHeel<>l"K>le*4ei|eH<HeHeHeHeHe 
%  i|ei|eHeHe*Hei|ee|ei|e*HeHei|ei|ce|eHee|e4ei|eeKi|eHe*i|e>l<el<e|e>|e*He4eHei|eHei|(*>|e 

%  HeHeHeHe>|eHes|eHeHeHei|ei|eHei|eHeHeHeHeHtHeHeHcHeHei|eHe]|ei|eHeHeHeHe>l'H<eKHeHcHeHe*HeHeHii|eHeHeHei|cHeHee|ee|eHeHe  10 

%  This  function  returns  the  value  for  a  ZD  integration 
%  for  the  slope  correlation  matrix  xy  and  yx  terms. 

%  The  integral  is  evaluated  for  several  values  of  y, 

%  then  sent  to  intgrate.m  to  do  the  ZD  integration. 

%  4ei|e4ei|ce|e>|e*HeHe>|eHe>|e>|ee|e>|e*ete**Hie|eHeHeHe**HeHe«*i|eef ****»*««*»****«««*«**« 

function  val  =  kolxyfun(xoflf,yoff) 

yliin=.5;  20 

step=.l; 

c=0; 

for  x=-ylim:step:ylim; 

c=c+li 

ma(c)=quad8('kolxy',yoff-ylim,yoiF+ylim,[],[],x,xoflf); 

end 

ma=ma'  i 
h=[step]; 

val=intgrate(ma,h);  30 


%  :HHttt********’lt*<S***'S****************** 

%  He*He*He*»*’|ee|«|ee|e*HeHei|BHci|e**«*3Kl|e***i|eHeHe4e*e|ee|‘*el‘* 

%  ♦*  Troy  B  Van  Caster  GE-97D  ** 

%  ♦*  Slope  Correlation  Matrix  for  ** 

%  **  xy  or  yx  directed  slopes  ** 

%  **  Last  Modified  5  Jun  97  ** 

%  +  +  ♦*♦*♦*♦*!(<  +  ***♦****♦***********♦♦*** 

%  HeHeHe>|eHe’lee|e*Hei|eHeHeHeHe«He*HeHei|eHeHe*He**HeHeHeHeHe****«i|e 

^  3|cHc3teHeite4e’teHeHeHcHeHee|eetee|eete4eettHeHe]teHeHeHcHeHeite3teHeHe])eHeHeH‘>teH(HcHe3te]te4:HeHeeKHeHeHeHcHcitei|eHeHeHc  XO 

%  This  function  is  called  by  kolxyfun.m  and  is  used 
%  to  evaluate  the  correlation  for  slopes  in  the  xy 
%  or  yx  direction  when  the  limits  of  integration  span 
%  the  origin.  The  modified  Bessel  function  uses  a 
%  two  term  expansion  to  prevent  the  integral  from 
%  having  an  infinite  value. 

%  HeHeHe>K>|ee|e>K*s|e*i|eH(H<HeH(i|e*i|eHeHe*eteH(I|eefeefei|eHes|e*i|ee|eHe**efe*i|eHeHe*HeHee|eHei|ei|eHei|ee|ee|eHei|ei|e 

function  val  =  kolxy(y,x,xoff) 
val=— 3.44*(. . . 

((y+l/2).-2+(-l/2+xoff-x)-2).-(5/6)  20 

((y+l/2).*2+(l/2+xoir-x)-2).-(5/6)  . . 

((y-l/2).'2+(-l/2+xoif-x)-2).-(5/6)  +•  • . 

((y-l/2).-2+(l/24-xoff-x)-2).-(5/6)  ); 
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E.2  Source  code  for  building  von  Karman  correlation  matrices 


function  [R]  =  vcorrmat(asize,d,ro,lo) 

%  He  3|e  H(  >((’((  if  lit  lit  4(>l«  Iff 

%  **  Troy  B  Van  Caster  GE-97D  ** 

%  *♦  von  Karman  slope  correlation  ** 

%  *♦  matrix  generator  ** 

%  **  Last  Modified  5  Jun  97  >i<h< 

%  steHtHtHt********************************* 

%  HtHtHcftHeHtHeiltHcHcHeHtHtHtHtHtHtHtHiHtHtHtHtHeHeittHtHtittittHiHtHtHtHtHcHt 
% 

%  This  builds  a  slope  covariance  matrix  using  the  von  Karman 
%  structure  function.  The  program  calculates  the  first  row 
%  values  then  uses  symmetry  to  fill  in  the  rest  of  the  matrix. 

%  The  program  builds  x  offset  and  y  offset  matricies  as  well 
%  as  a  ref  matrix,  xoff  and  yoff  are  subap*subap  in  size  and 
%  contain  normalized  x  and  y  offsets  from  xl,yl  to  x2,yS.  The 
%  ref  matrix  is  used  to  fill  in  the  rest  of  the  correlation 
%  matrix.  It  is  indexed  using  the  x  and  y  offsets  and  contains 
%  a  flag  to  the  appropriate  first  row  element  to  place  in  the 
%  matrix.  Called  by  R=vcorrmat(asize,d,ro,lo) 

%  This  program  calls  the  following  functions: 

% 

%  vonkxx  -  for  correlation  of  slopes  that  are  in  the  same  direction 
%  for  XX  slopes:  val=quad8(’vonkxx’,-l+xoff,l-hxoff, [],[], xoff, yoff) 

%  xoff  is  the  X  offset  normalized  by  d 

%  yoff  is  the  y  offset  normalized  by  d 

%  limits  of  integration  are  from  -1  to  1 

% 

%  vonkxyfun  -  returns  xy  slope  correlation 
%  calls  function  intgrate.m  and  vonkxy.m 

%  xoff  is  the  X  offset  normalized  by  d 

%  yoff  is  the  y  offset  normalized  by  d 

%  limits  of  integration  are  from  -112  to  112 

% 

%  vonkxy  -  determines  xy  slope  correlation 
%  val=quad8('vonkxy’,xoff-ll2,xoff+ll2,[  ],[ ],y, xoff, yoff) 

%  must  be  used  in  a  2D  integration  format 

%  HtHtHeHtHe>|t4:Ht«>|tHtH(Ht«it(it(H<ittH(HtHeHeHet1:HtH(He4(H(H(it(4<HtHtHtHeHeite9ftHeH(HtH(H(t|tHt)|(HeHtH(H(HtHcH(«HtHtHtH(HtHtHeHtHc>ttH( 

subap=asize'2;  %  length  of  covariance  matrix 

%  li*****************************************!**********^^!:^*********** 

%  build  X  and  y  offset  matricies  and  the  reference  matrix 

%  *♦*♦*+*♦*♦*♦*****♦******♦♦*********+♦**♦♦♦♦♦**♦**♦♦♦*♦****+*+***** 

index=0; 

for  x=l:asize 

for  y=l:aaize 

index=index+ 1 ; 

ref(y,x)=mdex*2— 1; 

xtemp(x,y)=(y-x); 

end 

end 

for  x=l:asize 
for  y=l:asize 

xoff(asize*x— asize+l:asize*x,asize*y— asize+l:asize*y)=xtemp; 

end 

endM 

for  x=l:asize'2 
for  y=l:asize‘2 

yoff(x,y)=(floor((y-l)/asize)-floor((x-l)/asize)); 

end 

end 

%  HtHtHiHt*Ht>|t*!|tHi>|tHt*HtHe>|tHtHtHiHiHtHtH:HtHtHiHtHi>|t***HtHt>|iHiHt*Htt|i*i|tHtHi*t|t*HtHtHiH<Ht**HtHiHi«HiHiHiHi>|t4tHi>|i 

H  calculate  the  first  row  of  the  correlation  matrix 
x=l; 

for  y=l:subap>i‘2 
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if  (floor(y/2))  "=  (y/2)  %it8  an  xx/yy  correlation 

R(x,y)=quad8('voiikxx',yoff(H-floor(x/2),l+floor(y/2))— l,yoff(l+fioor(x/2),l+floor(y/2))+l,. . 

[],[], xoff(l+floor(x/2),l+floor(y/2)),yoff(l+floor(x/2),l+floor(y/2)),lo,d)| 
else  %its  an  xy/yx  correlation 

%check  to  see  if  the  subapertrues  are  in  the  same  row/ column 
if  xoff(l+floor(x/2),floor(y/2))*yoif(l+floor(x/2),floor(y/2))==0 
R(x,y)=0; 
else 

R(x,y)=vonkxyfun(xofF(l+floor(x/2),floor(y/2)),yoff(l+lloor(x/2),floor(y/2)),lo,d); 

end 

end 

end 

%  Fill  in  the  rest  of  the  upper  diaganol  matrix 
%  start  with  y  rows 

for  x=2:2:subap*2 
for  y=x:subap*2 

if  (floor(y/2))  ==  (y/2)  %its  an  xx/yy  correlation 

R(x,y)=R(l,  ref(abs(yoff(floor(x/2),floor(y/2) ))+!,. . . 

abs(xoff(floor(x/2),floor(y/2)))+l)); 
else  %its  an  xy/yx  correlation 

R(x,y)=R(l,ref(abs(yoff(floor(x/2),floor(y/2)+l))+l,. . . 

abs(xoff(floor(x/2),floor(y/2)+l))-f-l  )+l); 

%check  the  x  offset  and  if  negative,  make  the  correlation  value  negative 
if  (  xoff(floor(x/2),floor(y/2)+l)  )*(  yoflF(floor(x/2),floor(y/2)+l)  )<  0 

R-(x.y)=-i*R(x.y); 

end 

end 

end 

end 

%  add  X  rows 
for  x=3:2:8ubap*2 
for  y=x:subap’t<2 

if  (floor(y/2))  '=  (y/2)  %its  an  xx/yy  correlation 

R(x,y)=R(l,  ref(ab8(xofF(floor(x/2)+l,floor(y/2)+l))+l,. .. 

ab8(yoff(floor(x/2)+l,floor(y/2)+l))+l  )); 
else  %its  an  xy/yx  correlation 

R(x,y)=R(l,ref(ab8(xoff(floor(x/2)+l,floor(y/2)))+l,. . . 

ab8(yoff(floor(x/2)+l,floor(y/2)))4-l)+l); 

%check  the  x  offset  and  if  negative,  make  the  correlation  value  negative 
if  (  xoff(floor(x/2)+l,floor(y/2))  )»(  yoiF(lIoor(x/2)+l,floor(y/2))  )  <  0 
R(x.y)=-i*R(x,y); 

end 

end 

end 

end 

%mirror  the  matrix 

R=R+triu(R,l) '  • 

R=R*(l/d)'2*(l/ro)-(5/3)| 


%  **  Troy  B  Van  Caster  GE-97D  ♦» 

%  **  Slope  Correlation  Matrix  for  ** 

%  **  XX  or  yy  directed  slopes  ** 

%  **  Last  Modified  3  May  97  ** 

% 

%  *i|<4c*!|c4iD!N»l<4i>l<*4<i|<4<***>l‘*t*4<*4<4<***4:!|i******>l<**4<*****’K*>l<4>**** 

%  This  function  is  called  by  vcorrmatm  and  is  used 
%  to  evaluate  the  correlation  for  slopes  in  the  xx 
%  or  yy  direction. 

%  The  true  function  is  evaluated  until  the  argument 


%  below  a  certain  level,  then  the  modified  Bessel 
%  function  uses  a  series  expansion.  There  is  also 
%  a  check  to  make  sure  that  the  series  expansion  is 
%  used  when  the  correlation  between  adjacent  subapertures. 

%  This  prevents  the  integral  from  taking  on  an  infinite  20 

%  value. 

% 

%  For  evaluation  in  the  xx  direction,  use  the  form: 

% 

%  quad8(’vonkxx’,yoff-l,yoff+l,[ ],[ ],yoff,xoff,lo,d) 

% 

%  For  evaluation  in  the  yy  direction,  use  the  form: 

% 

%  quad8(’vonkxx’,xoff-l,xoff+l,[ ],[],xoff,yoff,lo,d) 

function  val=vonkxx(y,xofiF,yoff,lo,d) 

%  evaluate  the  three  arguments 
argl=d/lo*(sqrt(xofF*2+y.‘2)); 
arg2=d/lo*(sqrt((xoff+l)‘2+y.‘2)); 
arg3=d/lo*(sqrt((xoff-l)‘2+y.*2)); 

threshhold=.02; 

40 

%  Define  some  initial  parameters 
A=0.08663»(lo)-(5/3); 

Bl=gamma(l/6)/pi*(l/6); 

B2=pi/(2*sin(5/6*pi)); 
trifun=(l-ab8(yofF-y))| 
terml=pl“(— 5/6)/gamma(l/6); 

%  Check  the  first  term  of  the  integral 

if  (argl  <=  threshhold)  |  (xoff<=l)  |(yoff<=l)  %  use  the  series  expansion 

term2=pi*(l+l/6)/gamma(l+l/6)*argl.*2;  50 

term3=pi*(5/6)/gamma(l+5/6)»argl. *(1+2/3); 
term4=pi*(2+5/6)/gamma(2+5/6)*argl. *(3+2/3); 
tl=trifun.*A.*Bl.*B2.*(terml+term2— term3-term4); 
else  %  use  the  actual  function 

tl=trifun.*A.»Bl.*argl.*(5/6).*besselk(5/6,2*pi*argl); 
end 

%  Check  the  second  term  of  the  integral 

if  (arg2  <=  threshhold)  |  (xoff<=l)  |(yofF<=l)  %  use  the  series  expansion 

term2=pi*(l+l/6)/gamma(l+l/6)*arg2.*2;  60 

term3=pi*(5/6)/gamma(l+5/6)*arg2. *(1+2/3); 
term4=pl*  (2+5/6) /gamma(2+5 /6)*arg2.  *  (3+2/3); 
t2=trifun.*A.*Bl.*B2.*(terml+term2— term3— term4); 
else  %  use  the  actual  function 

t2=trifun.*A.*Bl.*arg2.*(5/6).*bes8elk(5/6,2*pl*arg2); 
end 

%  Check  the  third  term  of  the  integral 

if  (arg3  <=  threshhold)  [  (xoff<=l)  |(yoff<=l)  %  use  the  series  expansion 

term2=pl*(l+l/6)/gamma(l+l/6)*arg3.*2;  70 

term3=pi*(5/6)/gamma(l+5/6)*arg3.  *(1+2/3); 
term4=pi*  (2+5/6) /gamma(2+5 /6)*arg3.  *  (3+2/3); 
t3=trifun.*A.*Bl.H<B2.*(terml+term2— term3— term4); 
else  %  use  the  actual  function 

t3=trifun.*A.*Bl.»arg3.*(5/6).*besseIk(5/6,2*pi*axg3); 
end 

val=2*tl— 12— 13; 


%  ******iitimi»*************************** 
%  NoK*********************************** 

%  **  Troy  B  Van  Caster  GE-97D  *♦ 
%  **  SD  integrator  for  xy  and  yx  ** 

%  **  slope  correlation  terms  ** 

%  *♦  Last  Modified  SO  April  97  ** 

%  ){(](( 3|c]|c :1c lie 

%  9lci)c3lcilc3lc!lc3lc]|c3lcc|c3|c9|cifc:lcilc3)c:lc](c9|c]tci|c)|c)|c:lc:lc:lci(c9|c])c))c]lcilc9tc)fc9tc))c))c 

% 
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%  This  function  returns  the  value  for  a  SD  integration 
%  of  the  slope  correlation  matrix  xy  and  yx  terms. 

%  The  integral  is  evaluated  for  several  values  of  y, 

%  then  sent  to  intgrate.m  to  do  the  SD  integration. 

function  val  =  vonkxyfun(xoff,yoff,lo,d); 
ylim=.5; 

step=.l;  20 

c=0; 

for  x=— yliin:step:ylim; 
c=c+l; 

ma(c)=quad8('vonkxy',yoff-ylim,yoflF+ylim,[],[],x,xofF,yoff,lo,d)i 

end 

ma=ma' ; 
h=[step]; 

val=intgrate(ma,h) ; 


E.3  Function  for  Two  dimensional  numerical  integration 


function  F  =  intgrate(f,h) 

%  function  F  —  intgrate(f,h) 

%  provides  numeric  integration  via  Newton  Coates  methods 
%  f  is  the  function  to  be  integrated  sampled  at  evenly  spaced 
%  intervals  h. 

%  f  -  Matrix  of  column  vectors  of  samples  of  functions  to  be  integrated 

%  h  ~  column  vector  of  the  width  of  the  spaces  for  each  column  of  f 

%  Let  N  =  length  of  sampled  function 
%  If  (N-l)l4  is  an  integer,  we  use  Newton  Coats  K=4 

%  If  (N-l)/S  is  an  integer,  we  use  Newton  Coats  K=S  10 

%  If  (N-lf/S  is  an  integer,  we  use  Newton  Coats  K=S  (Simpsons) 

%  If  (N-l)/l  is  an  integer,  we  use  trapezoidal  K=1  rule 

% 

%  K=l,  F  =  h/S*(f-l+f.N+sum(w*f(2:N-I));  w  =  [2  S  Z  . . .  2]; 

%  K=2,  F  =  h/3*(f.i+f-N+sum(w*f(2:N-l));  w  =  [4  2  4  2  ..  4] 

%  K=3,  F  =  3hl8*(f.l+f-N+sum(w*f(2:N-l));  w  =  [3  3  2  3  3  2  ..  3  3]; 

.  %  K=4,  F  =  2hl45*(7f.l+7f-N+sum(w*f(2:N-l));  w  =  {32  12  32  I4 

%  .  .32  12  32]; 

W  =  [1  2  1  0  0  0  0;  20 

1  4  2  1  0  0  0; 

1  3  3  2  1  0  0; 

7  32  12  32  14  7  0]; 

B  =  [1/2  1/3  3/8  2/45]; 

[N,M]  =  size(f); 

P  =  N-1; 
for  K  =  4:-l:l 

if  (abs(round(P/K)-P/K)  <  .0001) 

C  =  W(K,2:K+l)'*ones(l,round(P/K))  :  30 

b  =  B(K); 

A  =  reshape(C,ljN— 1); 

A  =  [W(K,1),  A(l:N-2),  W(K,K+2)]; 

F  =  (h*b)' : 
break; 
end; 
end; 
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